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Chapter  I 

nnxoDucriON 

With  the  introduction  and  proliferation  of  computers  into  all  facets 
of  the  work  place  and  the  home  environment,  a  new  awareness  of  the 
capabilities  and  short-comings  of  the  computer  for  various  tasks  has 
been  found.  The  computer  has  proven  very  useful  in  performing 
repetitive,  mundane  tasks  in  offices  and  manufacturing  process  control 
environments,  but  lack  of  a  good  real-world/computer  interface  prohibits 
many  uses.  Presently  a  computers  input  connections  to  the  real-world 
consist  mainly  of  a  keyboard  and  in  some  instances  joysticks,  graphics 
pads,  light  pens,  and  other  sensors  of  the  physical  world.  Recent 
research  into  this  interface  has  provided  the  computer  with  'ears',  that 
is  to  say  speech  recognition.  Not  only  can  the  computer  hear,  but  it 
can  also  act  upon  human  voice  commands  and  speech.  A  'voice'  and 
associated  language  generation  has  also  recently  become  a  reality.  The 
computer  can  generate  syntactically  correct  language  and  then  change 
this  into  intelligible  human  sounding  speech.  Perhaps  the  most 
isq>ortant,  and  by  far  the  most  complex,  interface  would  be  the  one  which 
gives  the  computer  'eyes'  or  sight.  Providing  the  computer  with  eyes 
and  vision  opens  new  realms  for  computer  automation  that  in  the  past 
were  either  too  difficult  to  perform  blindly  or  completely  impossible. 
This  vision  would  allow  the  computer  to  perform  difficult  and  tedious 
medical  or  industrial  inspections  at  a  much  higher  rate  than  is  possible 
with  human  workers  and  also  perform  the  inspections  under  hazardous  or 
harmful  enviromnental  circumstances.  With  knowledge  of  the  inspection, 
the  oomputer  could  also  make  decisions,  diagnoses  and  recommendations 
based  upon  its  previous  knowledge. 

The  possibilities  for  industrial  assembly,  combining  robots  with 
visual  feedback,  is  staggering.  With  such  devices,  the  jobs  thought  to 
be  too  boring  and  mundane  or  hazardous  in  the  manufacturing  field  can  be 
replaced  by  sighted  robots  that  do  not  complain  about  the  working 
conditions  or  monotony  of  the  job.  Just  as  computers  and  word 
processors  replaced  manual  record  keeping,  filing  and  typing  in  the 
modern  office  to  make  it  a  core  enjoyable  and  productive  working 


atmosphere,  the  seeing  computer  will  do  the  sene  for  the  industrial  and 
aannf ectnring  fields. 

A  few  of  the  many  possible  applications  for  intelligent  seeing 
computers  were  given  above,  but  by  no  means  even  scratched  the  surface 
as  far  as  the  nses  that  are  currently  envisioned  or  will  be  as  the 
research  evolves.  As  stated  above,  one  of  the  earliest,  and  cnrrently 
used,  applications  for  this  artificial  vision  is  in  the  area  of 
industrial  and  medical  examination  and  inspection 

[1] , [26] , [52] , [53] , [107] .  Same  of  the  inspections  are  those  which 
humans  may  deem  overly  tedious  and  boring,  while  other  inspection 
methods  performed  by  the  computer  can  not  be  carried  out  with  the  speed 
and  accuracy  required,  if  even  at  all,  by  human  workers.  On  a  scale  of 
difficulty,  that  of  the  visual  inspection  process  would  have  to  be  very 
near  the  lower  end  of  the  scale.  Very  little,  if  any ,  interaction  with 
the  article  being  inspected  takes  place.  Only  a  simple  keep/return 
deoision  or  perhaps  a  diagnosis  or  recommendation  may  be  required.  The 
problem  amounts  to  providing  the  system  with  a  'good'  product  or  example 
and  then  comparing  the  viewed  products  or  specimen  to  be  inspected 
against  the  good  model.  Provided  all  the  tests  and  requirements  are 
met,  the  piece  will  be  passed  or  a  good  diagnosis  returned.  If  any  one 
of  the  possible  visual  tests  fail,  the  article  can  then  be  rejected  and 
possibly  sent  back  for  correction.  In  the  medical  case  a  diagnosis  can 
be  generated  or  recommendations  for  further  tests  can  be  given. 

The  system  envisioned  here  is  good  for  its  particular  task  but  lacks 
the  intelligence  to  work  outside  the  narrow  area  for  which  it  was 
created.  Perhaps  more  simply  stated,  the  computer  can  not  'see* 
anything  that  it  has  not  been  programmed  to  see  and  hence  its  usefulness 
is  extremely  limited.  The  system  can  be  improved  with  provisions  to 
glean  information  from  its  field  of  view  that  represents  physical 
properties  of  the  objects  being  viewed  snd  their  relationship  to  the 
surrounding  area.  This  gained  intelligence  will  allow  the  computer  to 
make  judgements  about  its  environment  and  alter  it  to  meet  specific 
goals. 

This  would  allow  the  computer  to  perform  simple  manufacturing  and 
assembly  of  araltiple  parts  from  various  locations  snd  orientations. 
This  gain  in  ability  comes  at  a  high  cost  in  increased  complexity.  This 


system  is  aot  programmed  to  see  only  one  object  at  one  orientation,  but 
several  that  must  be  recognized,  tested  and  assembled  into  a  product. 
Even  this  problem  vould  have  to  be  classed  as  a  minor  one  when  compared 
with  the  visual  systea  that  humans  use  in  everyday  life. 

This  marriage  of  image  sensors  to  computers  has  already  started,  but 
the  hurdles  that  must  yet  be  crossed  to  reach  what  could  be  called 
'intelligent  machine  vision'  are  both  many  and  difficult.  Both  the 
computer  and  the  sensor  are  here  at  present,  but  like  a  young  child 
looking  at  the  world  he  has  never  seen  before,  can  not  understand  vhat 
he  sees  and  hence  cannot  constructively  affect  the  environment  in  which 
he  lives.  Ve  have  been  able  to  give  machines  the  eyes  and  brains  to 
see,  but  have  not  yet  been  able  to  give  them  intelligence,  understanding 
and  vision.  It  is  the  aim  of  this  work  to  provide  the  machine  with  a 
small  bit  of  information,  knowledge  and  understanding  so  that  one  day  we 
may  use  what  can  truly  be  termed  an  intelligent  seeing  machine. 


1.1  APPLICATIONS  FOR  ABTTFICIAL  VISION 

The  possibilities  and  end  uses  for  a  computer  that  can  intelligently 
see  are  endleas.  There  has  been  much  written  in  the  general  area  of 
machine  intelligence  and  more  specifically  machine  vision.  A  scan  of 
the  references  will  provide  many  articles  on  various  topics  dealing  with 
artificial  intelligence  and  vision.  At  present  the  systems  that 
actually  exist  and  are  in  use  are  few  and  far  between,  but  their  number 
is  growing  and  uses  expanding.  As  a  general  rule,  sighted  computers  can 
be  divided  into  two  major  categories.  The  first  category  consists  of 
systems  that  simply  inspect  their  respective  input  images  for  flaws  or 
abnormalities  [1] , [ 51] , [ 53] , [86] ,  [107],  and  [145]  and  make  decisions 

based  upon  these  investigations.  The  second  category  extends  the  first 
in  that  it  adds  feedback  to  change  its  environment  and  hence  act  upon 
what  it  has  seen  [61] ,  [67] ,  [77]  in  such  a  way  that  it  uses  vision  to 
constructively  interact  with  its  environment. 

Many  scientific  disciplines  have  sought  solutions  for  questions 
associated  with  vision.  Vhat  is  vision?  How  is  vision  defined?  How  is 
vision  interpreted?  These  are  just  a  few  of  the  endless  number  of 
questions  that  arise  when  dealing  with  the  topic  of  vision. 
Psychologists  look  to  the  internal  thought  processes  to  formulate 


theories  about  human  interpretation  of  vision  through  visual 
stimulation.  Physiologists  on  the  other  hand,  are  concerned  with  what 
biological  processes  are  necessary  for  the  generation  of  the  synaptic 
signals  associated  with  the  sensations  of  vision. 

From  research  carried  out  by  the  disciplines  above,  it  has  been  found 
that  the  human  visual  system  is  a  very  complex  and  complicated  network 
of  biological  subunits.  Seme  examples  of  these  subunits  are  the  light 
receptors,  i.e.  the  rods  and  cones,  the  optic  nerve  for  visual 
transmission,  and  the  brain  with  associated  memory  for  interpretation. 
Each  of  these  subunits  constitute  a  very  complex  system  in  and  of 
itself.  So  it  is  not  surprising  to  find  that  when  the  biological  visual 
system  is  modeled  by  hardware  devices  and  software  that  the 

nonbiological  visual  system  will  also  be  a  very  complex  and  complicated 
set  of  subunits. 

Just  as  the  human  visual  system  divides  the  whole  of  vision  into  many 
separate  but  related  parts,  so  then  should  an  artificial  vision  system. 
The  human  eye  is  not  simply  a  sensor  that  transforms  light  into 
electrical  signals,  but  also  an  enormous  parallel  processor  that 

performs  many  of  the  lower  level  visual  discrimination  functions.  Many 
of  the  rods  and  cones  of  the  retina  are  specially  sensitized  so  as  to 
fire  only  for  very  specific  physical  occurrences.  Examples  would  be 
those  that  fire  only  when  an  object  moves  at  a  specific  speed  in  a 
specific  direction,  those  that  fire  only  when  an  object  accelerates  and 
many  others  related  to  various  kinds  of  motion.  There  are  also  those 
which  fire  only  when  stimulated  by  an  input  field  containing  vertical  or 
horizontal  lines  or  boundaries.  With  this  it  can  be  seen  that  the  eye 
does  a  major  amount  of  processing  before  the  signal  even  enters  the 
optic  nerve.  The  same  cannot  be  said  about  the  camera  used  for  a 
computer's  eye.  By  the  time  the  visual  signal  reaches  the  brain  a  very 
large  amount  of  image  processing  has  already  taken  place.  The  brain 
receives  the  information  in  a  form  far  different  from  that  which  the  eye 
detected.  The  information  content  may  be  nearly  the  same,  but  the 
actual  amount  of  data  is  far  different.  The  sensors,  optic  nerve  and 
other  elements  that  connect  the  eye  to  the  brain  have  performed  what  is 
termed  data  compression. 
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1.2  IMPORTANCE  OP  IMAGE  CODING 

Of  mans  five  senses,  vision  vould  have  to  be  termed  the  most  data 
intensive.  Ye,  as  sighted  humans,  tend  to  take  the  amount  of 
information  processing  done  by  our  visual  system  for  granted.  As  was 
stated  in  the  previous  section,  the  amount  of  data  that  reaches  the 
brain  is  far  smaller  than  the  amount  of  data  that  is  actually  incident 
upon  the  retina.  This  data  reduction  or  compression  that  takes  place 
prevents  a  computational  overload  on  the  brain  by  delegating  some  of  the 
lower  level  visual  functions  to  other  parts  of  the  anatomy.  In  a 
similar  manner,  artificial  vision  needs  a  means  by  which  the  amount  of 
data  can  be  reduced.  It  is  at  this  point  that  the  similarity  in  the 
structure  of  the  human  and  artificial  visual  systems  must  stop.  The 
electronic  camera  is  a  much  simpler  device  than  is  the  human  eye,  and 
hence  processing  that  occurs  in  the  eye  must,  in  the  machine  model,  be 
centrally  processed.  For  this  reason,  and  perhaps  the  alleviation  of 
computational  complexity,  an  intelligent  visual  system  will  also  require 
a  means  for  data  compression. 

An  intelligent  vision  system  is  not  the  only  application  requiring 
the  need  for  data  reduction  of  imagery  data  and  computational 
requirements  are  not  the  only  resources  strained  by  the  very  large 
amount  of  data  present  in  imagery.  Storage,  transmission  and  system 
complexity  are  all  adversely  affected  by  the  high  data  rates  present. 
Many  ideas  remain  on  the  drawing  board  or  in  limited  use  because  of  the 
expense  that  this  large  data  rate  entails.  Hie  ideas  of  digital 
television,  video  phones  and  facsimile  have  been  with  us  for  sometime, 
but  the  enormous  bandwidths  associated  with  each  make  them  impractical 
when  it  is  realized  that  bandwidth,  like  many  material  resources,  is 
limited.  For  these  and  other  similar  products  it  can  be  seen  that  if 
commercial  viability  is  to  be  obtained,  then  the  video  bandwidth 
requirements  must  be  drastically  reduced  and  this  is  achieved  through 
data  compression. 

Applications  of  data  compression  are  not  limited  to  the  commercial 
market  alone.  The  military  and  space  scientists  also  have  many 
applications  that  for  one  reason  or  another  may  require  data 
compression.  On  earth  if  new  bandwidth  is  required  it  is  a  simple  task 
to  run  another  cable  or  optical  link,  but  for  space  it  is  much  more 


procedural  plots,  each  pair  consisting  of  two  different  gain  values. 
Figures  2.4,  2.6,  2.8,  and  2.10  consist  of  methods  which  require  forward 
and  inverse  transformations  at  each  iteration  while  the  remaining  plots 
depict  a  method  requiring  retransformation  only  on  block  intervals. 
Figures  2.4  through  2.7  involve  the  pel  recursive  displacement 
estimation  approach  while  2.8  through  2.11  involve  coefficient  recursive 
displacement  estimation.  All  of  the  plots  are  similarly  constructed 
where  the  horizontal,  or  I  axis,  represents  the  iteration  number  or  the 
number  of* times  equation  (2.12)  was  executed.  The  vertical,  or  Y  axis, 
repreruts  the  displacement  error  in  pixels.  The  test  images  used  were 
displaced  by  2  pixels  and  the  initial  guess  for  the  displacement  was 
assumed  to  be  zero,  hence  the  initial  displacement  error  of  2. 

The  data  compression  that  is  achieved  by  the  system  can  be  attributed 
to  two  major  characteristics  of  the  system  and  the  data.  The  first  is 
due  to  the  redundancy  removal  that  is  obtained  by  the  prediction  process 
and  the  second  is  dne  to  the  fact  that  many  of  the  transform 

coefficients  can  be  grossly  quantized  or  completely  neglected  with 

acceptable  results. 

One  of  the  major  problems  of  coefficient  recursive  estimation  is  that 
it  is  very  scene  and  position  dependent.  In  one  scene  or  position  the 

displacement  may  converge  in  as  few  as  4  or  5  iterations  where  for 

others  it  may  fail  to  converge  at  all.  Another  problem  may  prove  to  be 
the  choice  of  a  good  value  for  the  gain  factor,  epsilon.  As  figures  2.4 
through  2.11  verify,  the  choice  of  the  gain  factor  plays  an  important 
role  in  the  proper  convergence  of  the  algorithms  and  does  not  seem  to  be 
related  to  the  image  itself.  The  FORTRAN  computer  program  used  to  test 
the  algorithm  and  generate  the  iteration  plots  is  contained  in  appendix 
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Figure  2.3  Motion  Compensated  Hybrid  Transform-DPCM  coder-decoder  Pair. 


UN I  TAR 


Figure  2.2  Hybrid  Transform-DPCM  coder-decoder  Pair. 


Results  of  this  and  the  earlier  pel-recursive  method,  in  term  of  the 
displacement  error,  are  given  below  for  various  values  of  the  iteration 
gain.  The  test  image  is  a  radially  decaying  cosine  function  of  radius 
60  and  peak-to-peak  amplitude  of  220  at  the  center  decreasing  to  130  at 
the  circumference.  The  period  also  decreases  radially,  starting  with  a 
period  of  20  pixels  at  the  center  and  ending  with  10  pixels  at  the  edge. 
The  equation  used  to  generate  the  test  image  is 

I(R)  »  100exp(-0 .01R) cos(2nR/P)  +  128.  (2.14) 

The  displaced  image  is  shifted  two  pixels  in  the  x  direction  between 
time  frames.  Figure  2.1  is  a  picture  of  the  actual  test  image. 


FIGURE  2.1  Radially  Decaying  Cosine  Function 

Figures  2.2  and  2.3  point  out  the  differences  between  a  normal  hybrid 
transf orm-DPCM  coder-decoder  pair,  figure  2.2,  and  the  motion 
coaqiensated  hybrid  transf orm-DPCM  coder-decoder  pair,  given  in  figure 
2.3. 

Results  for  both  pel  recursive  and  coefficient  recursive  displacement 
estimation  for  a  single  block  from  the  test  image  are  contained  in 
figures  2.4  through  2.11.  The  plots  consist  of  4  pairs  of  similar 


These  transformed  blocks  are  then  arranged  into  a  col  nan  vector  by 
column  scanning  the  transformed  block.  They  define  the  n*h  coefficient 
of  the  block  of  the  transformed  image  to  be» 

ca(q)  -  IT(Xq.t)4a  (2.6) 

and  for  the  estimated  displaced  frame  from  the  previous  image. 

ca(q,B>  -  ITdq  -  D.t  -  t)*a.  (2.7) 

The  comparable  term  for  the  pel  recursive's  displaced  frame  difference 

A 

is  the  coefficient  prediction  error  ea(q.D)  and  is  given  by  the  equation 
below. 

•nU.D)  -  [I(Xq.t)  -  I(Xq  -  D.t  -  r)]Tia  (2.8) 

Minimization  with  respect  to  the  estimated  displacement  over  the  squared 
prediction  error  with  a  steepest  descent  form  is  given  below. 

Vl(«l)  -  V«l>  -<6/2)7Sa(q)ea(q.Da(q))  (2.9) 

Taking  the  required  derivative  will  yield  the  following  iteration 
formula. 

I>n+l(<l)  -  Dn  “  “#n(<1'°n(q))^Da{q),n(<l»°n(q))  (2-10) 

Rewriting  the  error  term  gives 

-®n-  «e(q-Da(<l))^a(q)(I(Vt)  “  I(Xq'Sa'  *"*>  >\J  (2n) 

A 

Noting  that  I(Xq, t)  is  not  a  function  of  DQ(q)  yields. 

Da+1(q)  ~  «e(q.Da(q))K7xIT(Xq-Da,t-T))in]  (2.12) 

When  going  to  the  next  block  the  initial  displacement  estimate  is  set  to 
the  final  estimate  of  the  previous  block. 

D0(q)  »  DN  N  _i(q-D  (2.13) 

r  c 

When  used  for  motion  compensated  interframe  hybrid  transform  DP  CM 
coding,  the  coder  transmits  a  quantized  version  of  the  coefficient 
prediction  error  whenever  it  exceeds  some  threshold.  This  allows 
updating  the  estimate  of  the  displacement  and  also  correcting  the 
prediction  coefficients. 
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This  aiy  tees  to  be  a  very  strict  requirement  in  tbit  very  little  reel 
world  notion  would  fit  the  model,  but  on  the  small  scale  of  the  pixel 
most  types  of  motion  can  be  approximated  by  nearly  pore  translation. 
The  background/ foreground  interaction  still  poses  convergence  problems. 
They  next  define  a  displaced  frame  difference  term  as  is  given  below  to 
represent  the  displacement  between  the  actual  value  and  its  estimated 
value. 


DFD<Xk,D)  -  I(Xk,t)  -  I(Xk  -  D.t  -  r)  (2.3) 

Where  in  this  case  the  term  S  is  defined  to  be  the  estimated  fora  of  the 
displacement  vector  D.  An  attempt  is  made  to  minimize  this  difference 
through  the  use  of  a  steepest  descent  algorithm  of  the  form  given  below. 


\+1  -  Dk  -  (e/2)7SklDFD(Xk.fik)]: 


(2.4) 


Here  e  is  a  gain  term  and  is  a  two-dimensional  gradient  operator 

a  h 

with  respect  to  Dk.  Taking  the  required  derivative  in  equation  (2.4) 
will  then  yield  the  following  iteration  formula. 


-  Dk  -  eDFD(Xk,Dk)pxI(Xk  -  Dk,t  -  t) 


(2.5) 


In  this  case  yz  is  the  two-dimensional  spatial  gradient  taken  with 
respect  to  the  row  and  column  directions  and  to  be  evaluated  at  the 

A  A 

point  X  ■  Xk  -  D.  This  term,  like  I(Xk  -  D.t  -  r)  .  may  require 

A 

interpolating  for  noninteger  values  of  D  at  each  new  iteration. 


When  actually  used  in  an  interframe  motion  compensated  image  coder, 

A 

the  transmitter  will  transmit  any  values  of  the  DFD(Xk,Dk)  term,  as  well 
as  the  required  address  information,  if  it  exceeds  some  set  threshold. 
This  quantized  correction  is  then  used  by  the  transmitter  and  receiver 
sections  to  update  the  appropriate  estimates. 


Following  this  same  motion  model,  Netravali  and  Stuller  [91] 
formulated  a  method  for  interframe  image  coding  termed  coefficient 
recursive  estimation.  It  is  an  extension  of  pel-recursive  with  the 
further  addition  of  a  unitary  transformation  so  that  the  operations  can 
be  carried  out  in  the  transform  domain.  The  two  methods  are  similar, 
assuming  the  same  motion  model,  but  now  the  image  is  partitioned  into 
rectangular  blocks  of  size  Nr  rows  by  NQ  columns.  Each  element  of  the 
block  is  then  multiplied  by  the  appropriate  transformation  vector. 
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Ohtsuka  [96]  define  a  similar  system  that  was  actually  implemented  in 
hardware.  The  method  nsed  for  determining  the  displacement  vector,  as 
before,  is  based  on  a  correlational  measure  but  now  also  takes  into 
account  the  displacement  vector  of  the  same  spatially  located  block  from 
the  previous  frame. 

Many  other  methods  have  been  employed  since  the  first  frame 
differencing  techniques.  A  method  by  Price.  Snyder,  and  Raj  ala  [111]  is 
based  upon  a  Fourier— Domain  filter  based  on  a  model  of  the  human  visual 
system  to  detect  motion.  The  model  breaks  up  visual  perception  into  two 
distinct  channels.  The  first,  their  so  called  z  channel,  is  a  temporal 
low  pass  and  spatial  band  pass  channel.  This  channel  carries  the 
information  contained  in  two-dimensional  patterns  with  high  spatial 

resolution  but  fairly  low  temporal  dependance.  They  believe  that  it  is 
this  channel  that  is  responsible  for  objects  with  structural  complexity 
but  with  little  or  no  motion.  The  other,  so  called  y  channel,  is  just 
the  opposite.  Its  characteristics  are  spatial  low  pass  and  temporal 
bandpass.  This  channel,  they  claim,  conveys  the  information  of  objects 
with  high  temporal  dependance  and  low  spatial  resolution.  It  is  this  y 
channel  that  is  used  for  the  detection  and  analysis  of  motion. 

2.2.3  Pel  Recursive  and  Coefficient  Recursive  Disnlacemen 

Stuller,  Netravali  and  Robbins  of  Bell  Laboratories  start  from  a 

completely  different  point  of  view  for  motion.  First,  the  end  product 
of  their  work  is  data  compression  and  not  tracking,  although  it  could  be 
modified  for  such.  The  model  is  used  for  normal  television  data  where 
adjacent  scan  rows  are  scanned  by  an  interleaved  method.  In  the  first 
method  of  pel-recursive  displacement  estimation  by  Robbins  and  Netravali 
[116],  the  image  model  for  pure  translation  with  no  background  is  given 
below. 

I(Xk,t)  =  I(Xk  -  D,t  -  r)  (2.2) 

Where:  I(Xk,t)  is  the  intensity  of  the  image  at  the  spatial  location  Xk 
and  time  t.  I(Xk  -  D,  t  -  x)  .  is  the  intensity  of  the  image  at  the 

spatial  location  Xk  adjusted  by  the  displacement  D  at  the  previous  time 

t-x .  This  image  model  is  defined  only  for  objects  undergoing  pure 
translation  and  not  involving  background  or  foreground  interactions. 
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be  possible  ss  wall  as  target  motion  prediction.  The  displacement 
vector  need  not  be  used  strictly  for  motion  related  studies  bat  may  also 
be  used  in  areas  of  data  compression,  remotely  piloted  vehicle  control, 
and  industrial  manufacturing  applications.  With  the  detection  and 
interpretation  of  stotion,  a  very  important  step  toward  artificial  vision 
will  have  been  obtained. 

Possible  solutions  of  this  problem  that  seem  nontractable  at  present 
may  in  the  near  future  be  made  posaible  due  to  the  advances  in  VLSI 
technology,  software  development,  parallel  processing  and  electro- 
optical  systems.  So  even  though  the  process  may  look  overly  complicated 
and  slow,  there  may  be  hope  in  the  future. 

2.2.2  Previons  Work  in  the  Field  of  Motion 

There  has  been  a  relatively  small  amount  of  work  carried  out  in  the 
area  of  machine  motion  analysis  until  very  recently  when  the  required 
hardware  and  software  became  available.  The  work  has  concentrated  in 
the  areas  of  displacement  estimation  and  interframe  image  coding.  One 
of  the  earliest,  and  perhaps  the  simplest,  method  used  for  motion 
detection  was  simple  image  frame  differencing.  That  is,  subtract  the 
previous  frame,  pixel  by  pixel,  from  the  current  frame  and  flag  as 
motion  any  difference  greater  than  some  set  threshold. 

M(i.j.t)  -  ABS[X( i, j , t)  -  I(i, j, t  -  t)]  (2.1) 

Motion  will  be  defined  to  have  occurred  whenever  M(i,j,t)  is  greater 
than  some  set  threshold.  Although  very  simple,  the  method  does  show 
good  results  for  a  very  limited  class  of  simple  images,  but  this  method 
has  many  drawbacks  that  will  limit  its  overall  usefulness.  First,  the 
output  is  very  sensitive  to  noise  in  the  input  images  because  it  is  a 
differentiating  type  process.  Also  any  camera  motion  between  image 
frames  will  translate  into  motion  at  every  pixel. 

Ihere  have  been  numerous  attempts  at  motion  compensated  image  coding 
in  the  past,  each  with  its  own  set  of  advantages  and  disadvantages.  An 
early  method  by  Oiorda  and  Raccin  [30]  breaks  the  image  up  into  a  number 
of  small  rectangular  blocks  and  determines  a  displacement  vector  for 
each  block  via  correlation.  If  the  displacement  can  be  found  it  is 
transmitted,  if  not  the  entire  block  must  be  transmitted.  Ninomiya  and 


previous  frame  and  if  possible  transmit  tbe  change  information  as  a 

function  of  motion.  Even  though  the  requirements  can  be  simply  stated 
the  actual  implementation  has  proved  to  be  difficult. 

2.2.1  Displacement  .gad  Motion 

Simply  stated.  motion  is  defined  to  be  a  time  series  of 

displacements.  That  is.  in  order  for  motion  to  be  perceived,  time  must 
pass  and  a  displacement  must  take  place.  If  artificial  viaion  and 
intelligence  is  ever  to  become  a  reality,  then  a  sufficiently  good  model 

for  motion  will  have  to  be  employed.  For  this  reason,  and  the  fact  that 

memory  space  is  always  limited,  the  vision  system  for  motion  analysis 
should  be  based  on  some  time  adaptive  displacement  algorithm. 

The  problem  can  now  be  stated:  Find  a  method  to  locate  the  spatial 
displacement  from  one  image  frame  to  the  next,  such  that  an  estimate  of 
the  direction  and  magnitude  of  any  detected  motion  can  be  made.  This 

estimate  should  be  based  upon  the  current  frame,  previous  frame  and 

previous  displacement  vectors. 

The  problem  statement  is  sifl^le  enough,  but  the  effort  is  complicated 
by  many  other  external  factors.  One  of  these  factors  involves 
object /background  and  obj ect/ foreground  interaction.  For  example,  if  an 
object  in  the  input  frame  moves  in  such  a  way  so  as  to  uncover  some  new 
background,  complications  will  arise  in  that  the  new  information  from 
frame  to  frame  now  consists  not  only  of  the  information  contained  in  the 

object  motion,  but  also  in  the  new  background  that  is  uncovered.  The 

imaging  system  should  have  the  capability  to  detect  the  difference 

between  the  moving  target  and  the  nomoving  uncovered  background. 

Another  problem  simpler  in  scope  than  the  above,  but  just  as  important, 
involves  the  loss  of  moving  objects  from  the  field  of  view  and  the 

addition  of  new  moving  objects  into  the  field  of  view.  Again  the  system 
should  be  able  to  detect  and  track  these  new  moving  objects  and  discard 
those  that  exit  the  field  of  view.  Finally,  objects  moving  in  the  field 
of  view  may  become  partially  or  totally  hidden  by  both  moving  and 

nonmoving  objects  and  the  system  should  be  capable  of  adapting  to  this. 

Once  the  motion  has  been  identified,  there  can  be  many  uses  for  the 
information  thus  provided.  Visual  tracking  of  moving  objects  will  then 


2.1.4  Other  Methods 


The  above  mentioned  methods  are  just  a  few  of  the  many  possibilities 
that  have  been  proposed.  Other  methods  include  normal  pulse  code 
modulation  (PCM)  as  well  as  adaptive  forms  of  PCM.  There  are  methods 
based  on  the  statistics  and  histogram  of  images  that  code  the  most 

common  gray  levels  into  short  code  words  and  rare  gray  levels  into  long 

code  words. 

Hybrid  coding  is  a  combination  of  both  transform  and  predictive 
coding  and  lately  has  been  studied  quite  extensively  [36].  [37],  [ 54 ]  . 
The  references  contain  many  other  variations  and  combinations  of  these 
and  many  more  methods. 

2.2  METHODS  EMPLOYING  MOTION  COMPENSATION 

Op  to  this  point  most  of  the  mentioned  coding  schemes  have  been  used 
for  intraframe  image  coding,  with  many  being  extended  to  include  image 
sequences  or  inter-frame  image  coding.  The  incluaion  of  the  third 

dimension,  temporal  or  time,  allowa  for  the  exploitation  of  the 
correlation  that  exists  in  this  temporal  direction. 

The  human  visual  system  is  a  very  fine  imaging  system  and  very 

susceptible  to  many  kinds  of  image  degradation.  The  degree  to  which 

certain  degradations  affect  the  imaging  process  vary  with  respect  to  the 
kind  of  degradation.  It  has  been  shown  that  the  eye  is  very  critical 
about  detail  in  still  pictures  but  becomes  less  critical  if  the  scene 
contains  motion.  If  the  entire  scene  is  undergoing  a  slow  constant 

translation  the  amount  of  detail  required  remains  high.  With  these 
simple  requirements  in  mind,  it  can  be  seen  that  any  inter-frame  coding 
system  with  a  human  as  the  final  viewer  will  be  required  to  maintain  a 
high  degree  of  detail  in  areas  of  little  or  no  motion  and  also  during 
times  when  the  camera  itself  moves  slowly  as  in  panning.  The  areas  that 
may  be  somewhat  degraded  by  the  coding  system  with  little  loss  of 

information  to  the  human  visual  system  are  the  areas  in  the  images 

undergoing  translation. 

The  area  of  motion  compensated  image  coding  lends  itself  directly  to 
these  requirements.  The  premise  of  motion  compensated  image  coding  is  a 
simple  one,  transmit  only  the  portions  of  the  image  that  differ  from  the 


the  predictive  techniques  can  not  offer.  Because  of  the  rather  high 
correlation  in  most  inages,  the  energy  tends  to  be  concentrated  in  the 
lower  frequency  or  sequency  coefficients.  Hence,  these  transfora 
coefficients  can  then  be  transnitted  if  their  energy  content  exceeds 
soae  value.  If  the  value  does  not  exceed  this  threshold  then  the 
coefficient  aay  be  grossly  quantized  or  not  transmitted  at  all.  Habibi 
and  Vintz  [35]  combined  this  method  with  block  quantization  to  obtain 
good  quality  results  on  the  order  of  2  bits/pixel  when  the  original 
image  was  coded  with  8,  or  a  data  reduction  by  a  factor  of  4.  Ngan  [95] 
extended  the  normal  transform  coding  with  the  addition  of  adaptation  to 
both  the  quantization  and  bit  selection. 

There  are  many  different  unitary  transforms  that  can  be  used  for 
transform  domain  image  coding.  Some  of  the  more  popular  are  the 
Hadamard,  Fourier,  Slant,  Cosine,  Sine  and  tarhunen-Loeve.  Each  method 
packs  the  energy  differently  and  different  results  will  be  obtained  if 
the  discard  method  remains  the  same.  Hideo  Kitaj ima  [62]  has  provided  a 
method  for  determining  the  energy  packing  efficiency  of  one  of  the  more 
useful  transforms,  the  Hadamard  transform. 

2.1.3  Predictive  Codin* 

Another  videly  used  coding  technique  for  digital  images  is  vhat  is 
generally  termed  predictive  coding  [19] , [23] , [36] , [37] , [38] , [55] , [92] , 
and  [122] .  The  predictive  coders  differ  from  the  above  two  methods  in 
that  data  is  not  just  discarded.  As  the  name  implies,  a  prediction  is 
made  by  using  local  pixel  intensity  values  in  some  predescribed 
combination  and  then  quantizing  and  transmitting  the  error.  Because  of 
the  high  degree  of  correlation  normally  occurring  between  adjacent 
pixels  the  prediction  process  works  quite  well,  which  implies  that  the 
error  is  normally  much  smaller  than  the  original  signal.  This  smaller 
error  signal,  in  terms  of  the  mean  square  power,  is  what  enables  the 
predictive  coder  to  obtain  its  high  degree  of  data  compression.  The 
more  structured  and  correlated  the  image,  the  better  the  predictor  is 
able  to  function  and  hence  the  lower  the  data  rate  for  a  fixed  fidelity 
criterion. 


2.1  PREVIOUS  METHODS  FOR  IMAGE  CODING 


Here  is  a  wealth  of  information  available  on  the  topic  of  general 
data  coding  and  picture  coding  in  the  literature.  Many  different 
methods,  applications  and  theories  are  available  from  a  large  number  of 
authors.  A  review  of  many  of  coding  solutions  is  available  from 

Ne travel i  and  Limb  [92].  Habibi  [35]  —  [38] ,  Pratt  [110].  Castleman  [14], 
Gonzales  and  Vintz  [31],  and  Andrews  and  Hunt  [4].  Some  of  the  more 
important  methods  will  be  briefly  described  in  the  following  sections. 

2.1.1  Sub semolina 

One  of  the  simplest  forms  of  data  compression  for  digital  images 

involves  subsampling  the  source  image  and  interpolation  of  the  discarded 
points  at  the  receiver.  This  can  be  accomplished  in  two  different  ways. 
Hie  first  involves  using  only  the  upper  right  hand  pixel  value  of  each 
of  the  N  by  N  sub-blocks  and  transmitting  that  value  to  represent  the 
entire  sub-block.  The  data  rate  here,  aa  in  the  second  method,  is  only 
l/N2  of  the  original  data  rate.  The  second  method  is  similar  to  the 
first  but  transmits  the  average  of  the  sub-block  instead  of  the  upper 
right  hand  pixel.  As  would  be  expected,  these  methods  tend  to  discard 
much  of  the  information  in  order  to  decrease  the  data  rate.  If  the 

image  was  already  sampled  above  the  Nyquist  rate,  then  some  data 
reduction  could  take  place  with  little  loss  to  the  quality  of  the 

output.  Normally  the  image  is  sampled  less  than  the  Nyquist  rate  and 
any  subsampling  will  cause  a  substantial  information  loss. 

2.1.2  Transform  Codlna 

The  simple  method  of  discarding  data  has  been  shown  not  to  function 
adequately  for  most  images,  but  with  the  inclusion  of  a  unitary 
transformation  data  can  be  discarded  with  much  less  loss  in  image 
quality.  The  function  of  the  unitary  transform  is  to  redistribute  the 
image  intensity  energy  in  such  a  way  that  it  is  mostly  contained  in  a 
small  number  of  transform  coefficients.  It  does  this  by  decorrelating 
the  spatial  data  and  hence  minimizing  the  statistical  redundancy  of  the 
information  to  be  coded.  Because  the  coding  is  done  in  parts  or 
sections,  an  error  in  one  part  will  not  propagate  to  other  sections  if, 
for  example,  a  channel  error  is  made.  This  is  an  advantage  that  many  of 


Chapter  II 

MOTION  COMPENSATED  TJtAGE  CODING  AND  DISPLACEMENT 
INTIMATION 

la  chapter  one  soae  of  the  requirements  and  restrictions  that 
necessitate  data  compression  as  well  as  some  possible  applications  for 
its  nse  were  presented.  In  this  chapter  the  basic  principles  of  data 
compression  will  be  reviewed  and  some  examples  of  various  methods  will 
be  presented.  Only  a  few  of  the  many  methods  proposed  in  the  past  will 
be  presented  dne  to  the  very  large  number  of  papers  written  on  the 
subj  ect. 

The  purpose  of  data  compression  is  to  take  a  data  sequence,  be  it 
image  data,  speech  or  any  other  information  source,  and  perform  a  data 
manipulation  in  such  a  way  so  as  to  be  able  to  reproduce  the  original 
signal  with  an  acceptable  amount  of  degradation  using  a  smaller 
bandwidth.  The  degree  to  which  this  goal  is  reached  is  based  both  upon 
the  statistics  and  form  of  the  image  as  well  as  the  operation  of  the 
coding  system  itself.  This  data  compression  is  often  schieved  through 
removal  of  the  large  amount  of  interpixel  correlation  that  is  normally 
present  in  most  data. 

One  of  the  early  uses  of  data  compression  was  in  the  field  of  digital 
voice  communication.  The  signal  generated  by  the  human  vocal  system 
tends  to  be  quite  correlated  and  hence  various  predictive  techniques 
have  been  used  to  obtain  good  compression  results.  When  the  field  of 
digital  image  coding  opened,  many  of  the  methods  used  for  speech  were 
directly  adapted  for  use  with  the  digital  images.  These  procedures  were 
not  able  to  take  full  advantage  of  the  geometrical,  statistical  or 
cognitive  structure  of  the  image  data.  Voice,  being  a  single 
dimensional  signal,  does  not  have  the  geometrical  structure  or  the  two- 
dimensional  statistical  dependence  that  image  data  has.  The  optimal 
compression  scheme  should  be  able  to  take  full  advantage  of  any 
structure  that  the  actual  image  has  to  offer.  In  most  cases  only  a 
small  subset  of  this  structure  is  actually  exploited. 


difficult  and  expensive  to  pat  up  a  new  satellite  or  add  the  required 
hardware  on  a  deep  space  probe.  Military  applications  can  include 
remotely  piloted  vehicles,  remote  surveillance  and  other  imagery 
applications.  One  side  effect  of  data  compression  that  often  times  is 
important  in  military  applications  is  that  of  secure  data  transmission. 
If  the  transmitter  and  receiver  are  not  a  matched  pair,  the  receiver 
will  not  be  able  to  recombine  the  data  sequence  to  obtain  the  original 
image.  There  are  many  different  ways  to  achieve  this  bandwidth 
reduction  and  most,  if  not  all,  exploit  the  high  correlation  between 
adjacent  pixels  to  reduce  the  data  rate. 

In  chapter  two  previous  work  in  the  field  of  image  coding  is 
presented.  It  provides  a  vide  overview  of  the  many  different  methods 
that  have  been  tested  and  employed  in  the  past  for  bandwidth  compression 
or  reduction.  The  emphasis  of  chapter  two  is  in  the  area  of  motion 
compensated  image  coding.  This  is  coding  that  exploits  knowledge  of  the 
motion  between  frames  in  sequential  image  data.  In  chapter  three  a  new 
method  for  motion  compensated  image  coding,  namely  Prediction 
Coefficient  Energy  Concentration  is  presented.  The  derivation  of  the 
displacement  estimation  procedure  as  well  as  the  coding  procedure  is 
presented.  In  chapter  four  information  on  the  implementation  and 
testing  of  the  Prediction  Coefficient  Energy  Concentration  model  is 
provided.  It  presents  the  information  that  ties  the  theoretical  aspects 
of  the  algorithm  to  its  actual  simulation  and  modelling.  The  results  of 
the  model  as  represented  by  a  software  implementation  are  presented  in 
chapter  five.  Actual  output  images  as  well  as  the  analysis  of  the 
algorithm  is  presented. 
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FIGURE  2.4  Modified  Pel  Recursive  Iteration  Graph 
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FIGURE  2.5  Pel  Recursive  Iteration  Graph 
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FIGURE  2.8  Modified  Coefficient  Recursive  Iteration  Graph 
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FIGURE  2.9  Coefficient  Recursive  Iteration  Graph 
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FIGURE  2.10  Modified  Coefficient  Recursive  Iteration  Graph 
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FIGURE  2.11  Coefficient  Recursive  Iteration  Graph 


2.2.4  Residual  Recursive 


The  method  of  Residual  Recursive  by  Rashid  and  Jones  [112],  [113]  is 
similar  to  both  the  pel  recursive  and  coefficient  recursive  methods  of 
the  last  section.  This  method  is  based  upon  an  earlier  procedure 
developed  by  Jones  [54],  [55]  for  image  coding,  namely  Adaptive  Hybrid 
Picture  Coding  or  AH  PC  for  short.  The  advantages  offered  by  the 
combination  of  AHPC  and  the  recursive  algorithms  developed  by  Stuller, 
Ne travail,  and  Robbins  lie  in  the  determination  of  the  gradient  term. 

Before  getting  to  the  motion  related  derivation,  it  is  necessary  to 
gain  some  background  information  on  AHPC.  AHPC  is  a  data  reduction 
coding  method  used  for  intraframe  data  compression.  The  method  is  based 
upon  a  two  step  correlation  removal  technique.  The  first  stage  of  the 
process  involves  a  one  dimensional  unitary  transformation  in  the  column 
direction  to  reduce  the  columnwise  correlation.  The  row  correlation  is 
then  reduced  through  a  prediction  process  that  operates  in  the  row 
direction.  Data  compression  is  achieved  because  the  mean  square  power 
of  the  prediction  error  and  predictor  overhead  is  less  than  the  mean 
square  power  of  the  original.  The  prediction  error  will  be  shown  to 
approximate  the  innovations  process  for  the  image  data  and  hence  the 
image  gradient.  A  block  diagram  for  the  AHPC  process  is  provided  in 
figure  2.12. 

The  transformed  pixels  are  modelled  as  a  one  dimensional 
autoregressive  series  where  the  predicted  value  of  the  signal  is  defined 
to  be  an  optimally  weighted  linear  combination  of  previously 
reconstructed  transformed  picture  elements  rQ.  The  weighting  factors 
are  the  predictor  coefficients  and  hence  the  prediction  equations  can  be 
written  as  follows. 


(2.15) 


The  difference  or  error  signal  eQ  is  then  defined  to  be  the  difference 
between  the  original  signal  sn  and  the  predicted  value  rQ. 


FIGURE  2.12  AnPC  Block  Diagram 


Using  the  *ean  square  error  as  the  optimality  criterion,  the  value 
for  the  predictor  coefficients  are  chosen  such  that  P#,  the  aean  square 
error,  is  ainiaized  and  given  belov. 


vJl.2 
•  **  _  1  ■ 


(2.17) 


Because  of  the  adaptability  of  the  algoritha,  the  predictor  coefficients 
a^  reaain  constant  for  a  single  row  or  learning  period  of  length  L.  Due 
to  the  statistical  change  of  the  data  fraa  row  to  row,  the  predictor 

coefficients  anst  also  change  at  row  boundaries.  Hence,  the  optiaal 

error  signal  for  a  single  learning  period  can  be  written  as  belov. 

e  -  s  -  AIa«  (2.18) 

The  saaple  value  n  is  allowed  to  vary  froa  1  to  the  learning  period  L. 
A*  is  a  col nan  vector  of  the  optiaal  set  of  predictor  coefficients  and 
Ka  is  a  row  vector  of  the  p  previous  reconstructed  transformed  picture 
eleaents. 

As  stated  earlier,  the  optiaal  predictor  vector  is  chosen  to  ainiaize 
the  average  error  signal  power  and  amounts  to  minimizing  the  following 
variance  tera. 

VAR(e)  -  kvAE(s)  -  COV(s,r)A  -  ATCOV(r,s)  +  ATVAR(r)A  ]  (2.19) 

L  o  n  &  & 

The  ainiaization  is  accoaplished  by  taking  the  partial  derivative  with 


respect  to  Aa. 


An(opt)  -  dVAR(e)/dAn 


(2.20) 


This  yields, 


Aa(opt)  -  VAR(r)’-lC0V(s) 


(2.21) 


When  complete  statistical  knowledge  of  the  source  is  not  known,  s  method 
is  needed  to  produce  good  estimates  of  the  statistics  during  each 
learning  period.  If  the  predictor  vector  is  treated  as  a  state  and  a 
linear  ainiaua  error  variance  sequential  estiaate  of  the  state  is 
desired,  then  sequential  estiaation  theory  can  be  used.  In  this  case 
Kalman  filtering  is  used. 


Given  the  optimal  set  of  predictor  coefficients,  the  reconstrncted 
transformed  picture  elements  can  be  generated  from  the  following 
equation. 

s_  *  A^R  (p)  +  e  (2.22) 

o  p  n  n 

Here  the  p  signifies  the  predictor  length.  RQ(p)  the  p  previous 
reconstrncted  transformed  picture  elements  and  is  the  quantized 

identified  predictor  vector  of  length  p.  The  identification  filter 
algorithm  from  [54]  can  then  be  written 

A(k+llS(k+l))  *  [I  -  K(k+l)RT(k) JA(klS(k) )  +  K(k+l)s(k).  (2.23) 

P 

The  difference  here  between  s(k)  and  S(k)  is  that  s(k)  is  a  single 
observation  while  S(k)  is  all  observations  up  to  and  including  s(k). 
Hence  the  desired  identification  algorithm  becomes. 

A(k+1)  =  A(k)  +  K(k+1) [s(k)  -  RT(k)A(k)l  (2.24) 

P 

where  the  bracketed  term  is  simply  the  prediction  error.  The  gain  term 
K(k+1)  is  given  by, 

K(k+1)  *  V. (k)R  (k)/(V  +  RT(k)V.(k)R  (k)).  (2.25) 

A  p  X  p  A  p 

VA(k)  is  siaqily  the  variance  of  the  identification  error  of  A(k)  as 
given  below. 

A 

A(k)  -  A(k)  -  A(k)  (2.26) 

Vz  is  a  scalar  term  not  obtainable  recursively  and  is  set  to  a  constant 
value.  Given  the  background  for  AHPC,  its  relationship  to  motion 
compensated  image  coding  can  be  discussed.  The  image  motion  model  for 
both  pel-recursive  and  coefficient  recursive  displacement  estimation  is 

I(Xk.t)  -  I(Xk  -  D.t  -  x)  (2.27) 

as  was  originally  given  in  equation  (2.2).  In  the  transform  domain  the 
equivalent  coefficients  are  given  below  as  originally  in  (2.6). 

cn(q)  -  IT(Xq,t)4n  (2.28) 

Note  the  values  given  for  the  displaced  frame  term  as  in  equation  (2.7). 
The  method  of  residual  recursive  displacement  estimation  lends  itself 
more  closely  to  coefficient  recursive  displacement  estimation  and  hence 


will  be  compared  as  such.  Note  first  that  the  transform  dona in 
coefficients  of  the  coefficient  recnrsive  method  are  also  the  signal 
samples  froa  AHPC. 


cn(q)  *  »„(q)  (2.29) 

That  is  to  say.  the  n**1  coefficient  of  the  q1*1  block  is  the  same  in  both 
notations  provided  the  linear  unitary  transform  is  the  same.  Hence,  the 
displaced  terms  then  become, 

°n(q.D)  -  sa(q,D)  (2.30) 

Replacing  this  notation  in  equation  (2.17)  yields 

•n(q)  -  «n(q)  -  A^6a(q).  (2.31) 

A 

Where,  as  before,  the  capital  letter  CQ(q)  defines  a  vector  of  the  p 
previous  predicted  values  froa  the  q^1  block.  On  a  coefficient-by¬ 
coefficient  basis,  this  can  be  written 


«n(q)  *  cn(<l)  "  kli*k®n-k(«)  * 


Likewise  the  displaced  term  becomes 


(2.32) 


,<q.D) 


o  (q,D)  - 


P 

.1 


VS)‘n-k(*'S) 


(2.33) 


Note  that  the  predictor  values  a^(D)  are  now  functions  also  of  the 

A 

estimated  displacement  D.  Recall  that  the  coefficient  prediction  error, 
now  with  the  prime  added,  is 


e^(q.D)  -  [I(Xq,t)  -  I(Xq  -  D,t  -  t)]Tia. 


(2.34) 


Which  can  also  be  written  as, 
e^(q,D)  -  oa(q)  -  ca(q,D). 


(2.35) 


The  comparable  term  in  residual  recursive  displacement  estimation  is  the 

A 

displaced  residual  difference  or  DRDQ(q,D) 


A  A 

DRDQ( q, D)  ■  ea(q)  -  ea(q,D). 


(2.36) 


As  before,  a  steepest  descent  algorithm  is  used  to  minimize  the  squared 

A 

difference  with  respect  to  DQ.  The  form  of  the  algorithm  is  given 
below. 


(2.37) 


Da+1  -  Dn  -(e/2)7g  [MU>  (q.B  )J2 
n 


Taking  the  indicated  gradient  results  in 


Vi  ■  D.  -  ■>*■>„ (?.■>„> 


(2.38) 


Before  determining  the  gradient  of  the  displaced  residual  difference 
ten,  it  is  necessary  to  rewrite  it  in  its  constituent  elemental  fon. 
From  (2.32), (2.33)  and  (2 .36) 

P  P 

D£Dn(q'V  "  cn(q)  "  J1Vn-k(q)  “  Ca(?,V  + Vk(q' V 

(2.39) 

The  cn_^(q)  i*  the  reconstructed  version  of  cn_^(q)  Noting  equation 
(2.3  5) 

P  P 

DRD  (q,D  )  -  e'(q.D  )  -  f  a.®  .(q)  +  f  a.  (D  )$  (q,D  )  (2.40) 

an  n  k  n-k  ^  k  n  n-k  n 

Rewriting  using  only  a  single  summation 

P 


DRDa(q.Dn)  -  o'(q.Da)  -  J^Vn-k^  "  W °n-k( q' V 


(2.41) 

At  this  point  it  is  necessary  to  make  use  of  the  following  identity. 

(2.42) 


n  n  n 

J  a  A  ■  JbB  - 
i-1  11  i-1  i-1 


That  is  given  a,  b.  A,  and  B  there  exists  a  set  of  c's  that  satisfy  the 
above  equation.  In  the  most  strict  sense,  if  they  pair  up  on  a  point  by 
point  basis  then  the  following  equation  will  hold. 


a^Aj  —  b^B^  —  c^(A^  “  B^) 


In  this  case. 


(2.43) 


ci  *  *  0>iBi)]/(Ai  -  Bt)  (2.44) 

Except  for  the  possible  case  where  A^  *  B^,  in  which  case  c^  will  equal 
zero.  So  with  this,  there  exists  a  set  of  coefficients  b^'s  that  will 
satisfy  the  following. 
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(2.45) 


DRDn(q.Dn) 


/4'V  “  ivV'Vk^1  -  ‘n-k^'V1 

k*l 


Because  of  the  good  prediction  and  error  quantization,  a  simplifying 
assumption  can  be  made  in  that  cfl  is  equal  to  cQ  or 

c  »  c  (2.46) 

n  n 

and 

cn(*'V  *  cn(q.Dn).  (2.47) 

With  this  equation,  (2.45)  becomes 


DRDa(q.Dn)  -  e^(q,D  )  -  f  bk(Da)  [ca_k( q)  -  cn_k(q.Da)]  (2.48) 

k-1 

Using  equation  (2.35)  yields. 


DRDn(<l.Dn) 


W*’ 0a)  -  ibk.^k(9.Dn) 
k«l 


(2.49) 


This  says  that  the  displaced  residual  term  is  composed  of  the  current 
residual  minus  a  linear  combination  of  previous  residuals  or  errors. 
Combining  (2.39)  with  (2.46)  and  (2.47)  yields. 


“Vl’V 


cn(q)  - 


c  <<1»D  >  - 
n  n 


P 

y  a  c 
kti  k  * 


Pa  a 

.  (q)  +  y  a  (D  )c  (q,D  ) 
i-k  k  n  n-k  n 

(2.50) 


A  A 

Taking  the  derivative  of  DRDn(q,Dfl)  with  respect  to  Dn  will  show  that 
the  first  and  third  terms  on  the  right  side  of  equation  (2.50)  will  go 

A 

to  zero.  These  are  not  functions  of  DQ.  The  remaining  gradient  term  is 
then. 


7A  DRD  (q,D  ) 

Dn  n  n 


-7a  c  ( q,D  ) 
Dan 
n 


C7A 

D 


n  k*l 


a  (D  )c  (q,D  ) 
k  n  n-k  n 


(2.51) 


This  is  the  most  general  form  of  the  gradient  equation,  note  that  the 
first  term  on  the  right  hand  side  of  (2.51)  is  what  is  termed  the 
coefficient  gradient  vector  in  coefficient  recursive  displacement 
estimation.  It  is  this  term  that  is  very  closely  approximated  by  the 
AHPC  residual  or  error  signal.  That  is. 
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(2.52) 


e  (q.D  )  *  c  (q.D  ) 
n  n  Dan 

a 

and  hence  it  is  not  required  to  calculate  the  coefficient  gradient  at 
every  iteration.  Replacing  this  in  equation  (2.51)  will  yield. 


5$  DRDn<q,Dn)  -  •  (q.$  >  £  ak(  V  °n-k(  V 

n  n  k=l 


(2.53) 


The  gradient  is  equal  to  the  AHPC  residual  plus  the  gradient  of  the 
prediction  tern.  Neglecting  the  second  term,  the  form  is  similar  to 
that  of  coefficient  recursive  displacement  estimation.  Using  equation 
(2.53)  in  equation  (2.38)  will  yield  the  following  iteration  formula. 

D_ . .  *  D  -  eDRD  (q,D  ) [e  (q.D  )  +  f  a.(D  )c  Jq.D  )]  (2.54) 

n+1  a  n  n  n  n  D  .A.  k  n  n-k  n 

n  k*l 

If  the  last  term  in  (2.54)  is  neglected,  the  following  iteration  formula 
is  very  close  to  the  form  used  in  equation  (2.12) . 


A  A  A  A 

D  . .  -  D  -  eDRD  (q.D  )e  (q.D  ) 
n+1  n  n  n.  n  n 


(2.55) 


Another  approach  that  used  all  of  the  available  information  .  i  be 
obtained  from  the  displaced  residual  difference  term  as  given  in 
equation  (2.36)  when  combined  with  (2.52)  to  yield  the  following. 


«V1- V  -  V»-'V 

n 

The  gradient  of  this  term  then  becomes, 

n  n 

and  the  iteration  formula  can  hence  be  written  is, 


(2.56) 


(2.57) 


(2.58) 


D-..  ■  ^  -  eDRD  (q,^  )V&  c  (q.D  ). 

a+i  n  n  nun  n 


(2.59) 


Results  for  this  method  are  given  in  figures  2.14  through  2.17  for 
the  same  image  data  as  was  used  for  the  earlier  methods.  The  FORTRAN 
program  used  to  test  the  algorithm  and  generate  the  plots  is  contained 
in  appendix  A. 
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FIGURE  2.13  Residual  Recursive  Displacement  Estimation 
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FIGURE  2.14  Residual  Recursive  Displacement  Estimation 
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FIGURE  2.16  Residual  Recursive  Displacement  Estimation 


function  of  th.0  data.  Hore.  as  in  tie  two  previous  methods,  the  minimum 
value  of  the  metric  or  cost  function  determines  the  nearest  integer 
displacement.  This  is  the  method  that  was  actually  employed  for  the 
results  presented  in  the  later  sections. 

3  .4  nTSPTJimiFNT  QUADRANT  ESTIMATION 

Given  an  estimate  for  the  integer  displacement,  the  next  problem  is 
then  to  determine  if  this  estimate  is  high  or  low  in  both  the  X  and  Y 
directions.  In  other  words,  the  displacement  needs  to  be  bounded  to  an 
integer  displacement  cell.  The  problem  can  be  better  visualized  by 
looking  at  the  four  displacement  cells,  or  quadrants,  about  the  integer 
estimate  as  is  shown  in  figure  3.5. 


FIGURE  3.5  Method  for  Minimum  Quadrant  Generation 

The  quadrants  are  numbered  I  through  IV  starting  with  the  upper  right 
hand  corner  and  moving  in  a  clockwise  direction.  Each  quadrant  has 
associated  with  it  four  values  of  the  metric,  one  being  the  minimum 
value  that  defined  the  integer  displacement  estimate.  Using  the  values 
generated  by  the  chosen  metric,  it  is  possible  to  determine  which  J 

quadrant  is  minimum.  This  is  accomplished  by  summing  the  values  of  the 
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Where  D  is  defined  to  be  an  integer  displacement,  1  is  the  actual 

A 

intensity  value  of  the  current  frame,  while  1(D)  is  the  current  frame 

A 

estimate  conditioned  on  the  integer  valued  displacement  vector  D. 
Therefore, 


R  - 


I  I  C(I(D),I)P.  A  .(1(D), I>). 
1(D)  I  11 


(3.29) 


The  joint  probability  on  I  and  1(D)  is  defined  as  ^(D) ,  j(I(D) ,  I)  . 
Rewriting  the  joint  density  in  terms  of  a  conditional  density  yields. 


R=  J  £  C(I(D).I)Pi|i(A)(l|l(D))Pi(A)(l(D)).  (3.30) 

All  of  the  terms  C(KD).I).  pi  |  i{D)  (I I H»> >  Pi(D)(I(°))  are  non“ 

negative,  so  R  can  be  minimized  by  minimizing  the  inner  sum. 

The  conditional  distribution  is  unknown.  It  is  therefore  assumed  to 
be  uniform.  This  is  equivalent  to  discounting  the  effects  of  the 
statistical  properties  of  the  image  sequence.  With  this  assumption,  it 
is  required  to  minimize 

R  »  £  C(I(D) , I) .  (3.31) 

I 


This  is  accomplished  by  determining  the  values  of  the  cost  function 
over  a  -p  to  p  neighborhood  in  both  spatial  directions  with  respect  to 
the  previous  frame.  Further,  a  two  valued  threshold  function  needs  to 
be  defined  for  the  cost  function. 

C(a,b)  *  0  if  b<a  (3.32) 

C(a.b)  *  1  if  b>a  f  (3.33) 

The  choice  of  this  metric  produces  a  uniform  cost  function.  When  this 
is  combined  with  the  absolute  value  method,  it  yields  a  counting 
arrangement  as  follows. 

N~1  N-l 

MU,j)  “IS  C(T,ABS[I(m,n.t)  -  I(aH-i,n+j,  t-r)  ] )  (3.34) 

“3)  n=0 

The  value  chosen  for  the  threshold  T  is  fairly  arbitrary,  but  somewhat 
loosely  related  to  the  variance  of  the  noise  and  may  even  itself  be  a 
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3 .3  .1  Isa  Si  Absolute  Values  si  Mqth<?4 


i 


The  first  aethod  is  based  on  the  sum  of  the  absolute  values  of  the 
difference  tera. 

N-l  N-l 

M(  i, j)  m  S  £  ABS[I(m,n,t)  -  I(m+i, n+j, t-r) ]  (3.26) 

a-0  n-0 

Where  both  i  and  j  are  allowed  to  vary  from  -p  to  p.  ABS  in  the  above 
iaplies  taking  of  the  absolute  value  of  the  quantity  in  the  square 
brackets.  The  best  integer  displacement  is  then  defined  to  occur  at  the 
point  at  which  the  aetric  M(i,j)  is  ainimua.  From  a  cost  function  point 
of  view,  the  absolute  value  function  defines  an  absolute  error  cost 
function. 

3 .3 .2  Isa  2l  Squares  si  Difference  Method 

The  aethod  is  similarly  defined  as  that  in  (3.26)  with  the  use  of  a 
squaring  function  instead  of  the  absolute  value.  This  offers  the 
advantage  of  penalizing  the  aetric  greater  for  large  differences  than 
for  a  number  of  saaller  differences. 

N-l  N-l 

M(i.j)  -  I  I  [I(a.n.t)  -  I(*+i,  n+j ,  t-c)  ]2  (3.27) 

s*-0  n-0 

Here  again  the  ainiaun  value  of  the  aetric  defines  the  location  of  the 
best  integer  estimate.  From  a  cost  funotion  point  of  view,  the  squaring 
function  defines  a  square-error  cost  function. 

3.3.3  Threshold  Count lnz  Method 

To  this  point  the  aethod  of  threshold  counting  has  proved  to  be  the 
aost  advantageous.  It  offers  the  advantage  of  not  penalizing  the  metric 
for  saall  differences,  on  the  order  of  the  noise  that  may  be  present, 
and  counting  all  values  above  this  threshold  equally.  As  in  the  two 
methods  above,  either  the  absolute  value  or  squared  difference  method 
could  be  used,  but  due  to  its  simplicity  the  absolute  value  method  has 
been  chosen.  The  procedure  utilized  produces  an  estimate  that  minimizes 
the  expected  value  of  the  cost.  The  expected  value  of  the  cost  is 
defined  as  the  risk.  That  is, 

S  -  E[C(I(D) ,1) ) ]  (3.28) 


'i 
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3.3 


The  method  outlined  above  requires  s  good  estimate  of  the  integer 
displacement  bounds.  Several  different  methods  have  been  tried  for  the 
purpose  of  establishing  the  displacement  and  a  few  of  the  various 
functional  metrics,  or  cost  functions,  that  can  be  used  will  be 
discussed  below. 

In  each  of  the  metrics  described  below,  the  previous  frame  search 
area  is  defined  over  a  -p  to  p  neighborhood  for  both  spatial  directions. 
The  metric  differencing  is  performed  on  a  pixel  by  pixel  basi  every 

pixel  in  the  current  block.  A  new  value  of  the  metric  is  generated 
every  time  the  current  block  is  shifted  with  respect  to  the  previous 
frame  block.  The  current  block  of  size  NxN  is  compared  to  all 
contiguous  NxN  subsets  of  the  previous  frame  local  neighborhood  of  size 
(2p+N)  by  (2p+N)  as  is  shown  in  figure  3.4.  Note  that  the  metric 
generated  will  be  a  (2p+l)  by  (2p+l)  matrix,  where  each  entry's  location 
is  a  function  of  displacement. 


The  problem  of  finding  the  noninteger  portion  of  the  displacement  has 
been  replaced  by  a  two  step  process  that  is  able  to  perform  a  similar 
function.  First,  for  the  one  dimensional  example,  the  actual 
displacement  must  be  bounded  above  and  below  by  a  high  and  low  integer 
estimate.  Secondly  the  two  values  for  the  delta  functions  must  be 
generated  via  the  solution  of  the  regression  equation  given  below. 

1 

f(nT)  -  £  a(i)  g(nT  -  L  +  i)  +  e(nT) .  (3.22) 

i«0 

Where  L  is  the  lower  integer  bound  of  the  displacement  estimate.  Given 
a  sufficiently  large  number  of  values  for  f(nT)  and  g(nT)  the  values  for 
a(0)  and  a(l)  can  be  found  through  least  squares  estimation.  Xhe  same 
procedure  can  be  extended  to  the  two-dimensional  process,  where  four  on- 
sample  delta  functions  replace  the  single  off-sample  function.  As 
before  note 

Ki.j.t)  -  6(x,y)*I(i,  j.t-r)  .  (3.23) 

where  x  and  y  denote  the  location  of  the  delta  function  which  need  not 
fall  on  the  two-dimensional  sampling  grid.  A  time-modified 
autoregressive  model  taking  into  account  that  the  prediction  is  made 
over  the  time  boundary  is  then  given  by 

P  P 

Ki.j.t)  *  's  y  a(m,n)I(  i-m,  j-n,  t— c)  +  e(i,j,t).  (3.24) 

■— P  n— P 

If,  as  before,  the  best  integer  estimate  of  the  displacement  bounds  can 
be  found,  a  displacement-updated,  time-modified  autoregressive  model  can 
be  given  as 

1  1 

Ki.j.t)  ■  £  £  a(m,n)I(  i-K+m,  j-L+n,  t-r)  +  e(i,j,t).  (3.25) 

m*0  n“0 

Here  K  and  L  designate  the  lower  bounds  of  the  integer  displacement 
estimate.  Here,  as  in  the  one  dimensional  case,  the  regression 
coefficients  are  able  to  perform  the  same  function  as  the  noninteger 
portion  of  the  displacement  vector. 
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For  example,  if  the  true  rain*  of  the  delta  function  occurs  at  t  * 
1.3  5T,  it  caa  be  replaced  by  two  differeat  delta  faactioaa.  One 

positioaed  at  t  -  IT.  aad  the  other  at  t  *  2T.  The  relative  value  of 

each  delta  fuactioa  is  obtaiaed  from  the  following  formulas,  where  &^(t) 

is  a  lower  bound  of  the  shift  fuactioa  aad  ^(t)  is  an  upper  bound  of 

the  true  value  of  the  shift. 

5j(t)  -  1  -  e  (3.20) 

S^t)  -  e  (3.21) 

la  the  example  given  above  the  value  of  the  fuactioa  at  t~lT.  labelled 

&}(t).  would  be  .65,  while  the  value  at  t*2T,  labelled  6^(0,  would  be 

.35  and  heace  is  able  to  perform  the  fuactioa  of  linear  interpolation  as 
if  it  were  a  noninteger  shift.  Figure  3.3(b)  shows  the  placement  of  the 
two  off  sample  delta  functions  for  the  aoaiateger  shift  ia  3.3(a). 


FIGURE  3.3(b)  Delta  Function  for  Noninteger  Shift 

The  problem  at  hand  uses  this  shifting  and  interpolation  ability  of 
multiple  delta  functions  to  bypass  the  need  for  sub-pixel  determination 
of  the  shift  vector.  Assuming  that  a  method  can  be  found  to  determine 
the  nearest  integer  displacement,  with  the  use  of  regression  analysis 
there  is  no  need  to  find  the  noninteger  portion  of  the  displacement. 
Instead,  only  the  magnitudes  and  locations  of  the  delta  functions  are 
required.  How  the  integer  displacement  is  determined  is  given  in  the 
next  section. 
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the  use  of  linear  interpolation  of  the  two  adjacent  points.  The  fora  of 
the  linear  interpolation  is  given  below  where  f^(n'T)  designates  the 
function  value  at  the  lower  bound  of  the  interpolation,  fu((i'+l)T) 
designates  the  function  value  of  the  upper  bound,  and  fm(n'T)  is  the 
function  value  at  the  linear  interpolated  point  between  these  bounds. 
The  n*  is  used  here  to  denote  a  fixed  value  for  n. 

fa(n'T+e)  -  fb(n'T)  +  e(ftt((n'+l)T)  -  fb(n'T))  (3.16) 

In  the  above  equation  e  is  the  distance  from  n'T  to  the  desired 

interpolation  point,  or  if  yon  wish,  the  noninteger  portion  of  the 

displacement.  Carrying  this  idea  further,  an  interpolated  function  can 
be  generated  by  letting  n*  vary  over  the  appropriate  limits  to  include 
all  values  of  the  sampled  function.  Rewriting  (3.16)  in  function  form 
yielda, 

fa(nT+e)  -  (l-e)fb(nT)  +  ef n( <n+l)T) .  (3.17) 

Note  that  in  the  two  equations  above,  the  shift  or  interpolation  is 
limited  to  the  range  of  0  to  1.  That  is,  e  is  bounded  to  the  interval  0 
to  1.  From  equations  (3.16)  and  (3.17),  it  can  be  seen  that  the 

interpolated  portion  of  the  function  is  a  linear  combination  of  the 

points  on  either  side.  It  is  from  this  equation  that  the  idea  of 
sraltiple  delta  functions  performing  interpolation  can  be  found.  The 
first  term  on  the  right  side  of  equation  (3.17)  can  be  modified  by  a 
delta  function  at  the  origin  with  value  1-e.  The  second  part  of  (3.17) 
is  modified  by  another  delta  function  but  at  location  tmT  and  of  value 
e.  Rewriting  (3.17) in  a  convolutional  form  yields, 

f_(nT+e)  -  ( 1-e) f(nT)*5( t)  +  ef  (nD*6(t-T)  .  (3.18) 

■ 

Hence,  the  single  off-sample  delta  function  located  at  t^e  has  been 
replaced  by  two  on-sample  delta  functions.  This  can  be  further  extended 
to  perform  for  any  value  of  interpolation.  That  is,  the  displacement 
does  not  have  to  be  restricted  to  the  0-1  range.  When  this  is  the  case, 
e  remains  the  noninteger  portion  of  the  displacement  and  both  delta 
functions  are  then  shifted  by  the  remaining  integer  amount  of  the 
interpolation  or  shift.  With  S  being  the  integer  amount  of  shift, 
(3.18)  can  be  written  as  follows. 


f-(nT+e)  -  (1-e) f(nT)*6(t-ST)  +  ef (nT)*6(t-(S-l)T) 


(3.19) 


a 

I 


FIGURE  3.2(b)  Delta  Function  fot  Integer  Shift 

function  can  be  obtained  by  a  linear  coabination  of  delta  functionc 
bordering  the  trne  time  or  phase  shift.  Dsing  the  method  of  linear 
interpolation,  the  single  off-sample  delta  function  can  be  replaced  by 
two  omrsample  delta  functions.  The  magnitude  of  each  delta  function  is 
based  on  the  location  of  the  actual  off-sample  delta  funotion. 

aCnT3  fCnT5 


FIGURE  3.3(a)  Sampled  Function  Noninteger  Shift 

Before  the  shift  can  be  determined  it  is  necessary  to  explain  the 
mechanics  of  this  linear  interpolation  process.  It  is  assumed  that  any 
value  of  the  function  between  sample  points  can  be  determined  through 


Hie  graphical  representation  of  this  case  is  given  in  fignre  3.2(a)  and 
hence  f(nT)  and  g(nT)  for  this  specific  exaaple  are  given  below. 


f(nT)  «  AcosCvnT  +  9)  (3.13) 

g(nT)  *  Acoa(w(nT  -  x)  +  9)  (3.14) 


FIGORE  3.2(a)  Saapled  Function  Integer  Shift 

As  long  as  x  is  an  integer  aultiple  of  the  saapling  period  T,  as 
shown  in. fignre  3.2(b).  the  same  holds  true  for  the  sasipled  case  as  for 
the  continuous  case  and  hence. 

f (nT)  -  g(nT)*8( t  -  are)  (3.15) 

where  arc  in  an  integer  aultiple  of  T. 

Figure  3.3(a)  graphically  depicts  the  usual  case  in  that  the  shift  is 
a  noninteger  aultiple  of  the  saapling  period.  The  problem  here  is  that 
&(t~r)  aay  not  fall  on  a  saapled  time  interval.  Exaaple  t  •  1.5T  or  1.5 
tiaes  the  saapling  interval  and  hence  would  lie  halfway  between  the  two 
saapled  tiae  intervals.  The  method  presented  here  attacks  this  problem 
of  noninteger  shift  with  the  help  of  linear  interpolation.  Assuming 
that  all  points  of  the  saapled  time  function  between  samples  can  be 
determined  sufficiently  close  by  linear  interpolation  of  the  sample 
points  on  each  side,  then  a  sufficiently  good  estimate  of  the  shifted 


g(t)  ■  Acos (»( t  -  r)  +  9) 

Equation  (3.4)  for  this  specific  example  then  becomes, 
f(t)  «  Acos (»( t  -  x)  +  9)*8(t— c) 

There  8(t— c)  is  shown  in  figure  3.1(b). 
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FIGURE  3.1(b)  Delta  Function  for  Continuous  Shift 

Using  the  Fourier  transform  to  perfora  tine  convolution  yields 

f(t)  -  F'1[FIAcos(«(t  -  r)  +  G)]  F[5(t  -  r)]]  (3.8) 

f(t)  -  Acos (u»t  +  9)  (3.9) 

where  F  designates  the  forward  Fourier  transform  and  f-1  the  inverse 
transform.  This  is  the  case  as  was  given  in  equation  (3.5) . 

Getting  closer  to  the  problem  at  hand,  the  sampled  versions  of  f(t) 
and  g(t).  namely  f(nT)  and  g(nT).  obtained  from  the  sampling  function, 

f (nT)  -  f ( t) y  8( t-nT)  (3.10) 

can  be  written  as  follows. 

g(nT)  -  f(nT-c)  (3.11) 


or 


f (nT)  -  g(aT)*6(nT) . 


(3.12) 


examining  a  single  dimensional  continnona  time  function  example  with  the 
help  of  linear  systems  theory.  The  two-dimensional,  noiseless,  no 

background  case  becomes,  in  a  single  dimension,  the  problem  of 

determining  the  phase  shift  between  two  identical,  but  time  shifted, 
time  functions.  For  example  if 

g(t)  -  f(t-r)  (3.3) 

the  two  functions  are  said  to  be  identical,  with  the  exception  of  the 
time  shift.  The  problem  then  of  finding  the  time  shift  amounts  to 

finding  the  value  for  x.  Using  the  concept  of  convolution,  g(t)  and 
f(t— c)  can  be  related  further  by 

f(t)  -  g(t)*5(t-t)  (3.4) 

where  the  function  *  is  defined  to  be  convolution  and  6(t-x)  is  a  delta 
function  located  at  t*t .  The  problem  remains  to  determine  the  time 

displacement  from  the  origin  to  the  location  of  the  delta  function. 
Figure  3.1(a)  graphically  preaents  a  possible  representation  of  these 
two  continnona  time  functions  shifted  in  time  by  r.  This  displacement 
can  be  determined  through  the  use  of  correlation  or  a  matched  filter 
such  that  the  output  is  the  phase  thift  between  the  two  signals. 


In  the  aethods  previously  discussed  for  notion  compensated  image 
coding.  [91].  [112],  [113],  snd  [116],  the  current  image  pixel  is 

modelled  as  a  pixel  from  the  previous  frame  displaced  by  some 
displacement  vector  D.  This  can  be  shown  as  in  (3.1). 

I(Xk,t)  -  I(Xk  -  D,  t  -  t)  (3.1) 

Vhere  I(Xk,  t)  is  the  pixel  intensity  at  location  Xk  and  x  is  the  time 
delay  between  adjacent  frames.  This  functions  adequately  provided  the 
displacement  is  of  integer  order.  That  is,  no  partial  pixel 
displacements  are  allowed.  If  sub-pixel  displacements  do  exist,  this 
necessitates  the  finding  of  a  D  vector  that  is  also  of  sub-pixel  order. 
A  somewhat  more  complicated  approach,  but  offering  the  advantage  of  not 
having  to  determine  a  displacement  vector,  can  be  found.  The  current 
frame  pixel  intensity  is  defined  to  be  a  linear  combination  of  pixel 
intensities  from  the  previous  frame. 

M  N 

I(Xk.t)  ■  '£  £  a(m,n)I(Xk(m,n) ,t— r)  +  e(Xk.t)  (3.2) 

m-1  n-1 

The  disadvantage  here  is  that  although  no  motion  vector  is  required,  a 
set  of  predictor  coefficients  a(m, n)  is  required.  If  the  system  is  to 
allow  for  a  p  pixel  shift  in  any  direction,  the  size  of  this  'a'  matrix 
would  be  at  least  (2p+l)  by  (2p+l).  Hence,  the  amount  of  data 
compression  that  can  be  achieved  is  greatly  diminished  as  p  gets  large. 
Note  that  this  prediction  matrix  must  be  transmitted  every  time  the 
displacement  changes  between  blocks. 

The  ideal  situation  would  be  to  take  advantage  of  the  linear 
combination  method,  so  that  the  accuracy  of  the  motion  vector  can  be 
kept  to  integer  displacements,  and  yet  exploit  the  data  compression  that 
the  displacement  vector  approach  offers.  Going  back  to  (3.2)  it  can  be 
seen  that  the  solution  of  the  problem  involves  the  determination  of  the 
prediction  matrix  a(m,n).  This  results  in  a  two-dimensional  regression 
problem.  Noting  what  the  physical  implication-  are  in  relation  to  the 
regression  problem,  it  can  be  seen  that  the  regression  coefficients 
would  be  used  to  perform  translation  and  hence  the  model  can  be 
simplified.  A  better  understanding  of  this  concept  can  be  obtained  by 


Before  getting  into  the  actual  algorithm,  there  are  a  set  of 
assumptions  and  requirements  that  need  to  be  stated. 

1.  Each  image  in  the  sequence  is  broken  up  into  a  set  of  square 
blocks.  The  motion  for  that  block  is  assumed  to  be  constant  over 
the  entire  block,  regardless  of  the  chosen  blocksize. 

2.  The  motion  as  modelled  by  the  algorithm  is  pure  translation,  that 
is  only  motion  that  is  translation  or  can  be  modelled  by 
translation  between  frames  is  actually  modelled.  Nontranslation 
type  motion  and  problems  of  occlusion  will  be  handled  by  the 
quantization  of  the  residual  and  not  entirely  by  the  displacement 
vector.  Although  in  some  cases,  the  prediction  coefficients  can 
somewhat  compensate  for  nonideal  motion. 

3.  It  is  further  assumed  that  the  maximum  displacement  between 
frames  is  known.  The  system  is  no  more  complex  for  large 
allowable  displacements,  but  the  number  of  calculations  increases 
rapidly  as  the  maximum  allowable  displacement  increases.  Hence, 
the  maximum  displacement  estimate  should  be  kept  as  small  as 
possible  to  improve  performance.  For  a  majority  of  image  work  a 
value  of  five  or  six  pixels  is  appropriate. 

4.  It  is  known  that  the  human  visual  system  will  accept  a  higher 

degree  of  degradation  for  scenes  undergoing  translation  than  it 
will  for  static  scenes.  For  this  reason,  the  blocks  that  are 
undergoing  translation  are  allowed  a  lower  signal  to  noise  ratio 
than  those  which  do  not  move.  That  is,  the  moving  blocks  can 

have  a  higher  degree  of  degradation  than  the  stationary  ones. 

With  this  set  of  assumptions  in  hand,  the  goals  of  the  algorithm  can 
be  stated:  the  overall  system  should  be  able  to  obtain  a  fairly  high 
degree  of  data  compression,  but  not  at  the  expense  of  overall  picture 
quality.  The  system  should  be  able  to  achieve  a  very  low  data  rate  for 
image  sequences  with  little  or  no  motion.  For  the  blocks  in  the  image 
that  are  undergoing  translation  the  method  should  be  able  to  determine 
the  motion  vector  to  sub-pixel  accuracy  and  code  the  resulting  output  at 
a  sufficiently  low  data  rate. 


Chapter  III 

PREDICTION  COEFFICIENT  ENERGY  CONCENTRATION 

In  the  last  chapter  some  of  the  previous  attempts  at  motion 
estimation,  and  in  particular  motion  compensated  image  coding,  were 
presented,  setting  the  groundwork  material  for  this  chapter.  The  method 
of  'Prediction  Coefficient  Energy  Concentration'  differs  from  the 
previous  methods  not  in  terms  of  optimal  output,  but  only  in  the  way  in 
which  this  goal  is  achieved. 

As  the  name  'Motion  Compensated  Image  Coding*  implies,  motion  or 
movement  estimation  is  used  to  decrease  the  data  rate  required  in  the 
transmission  of  time  sequential  image  data.  The  use  of  this  motion  or 
movement  requires  that  a  good  estimate  for  the  displacement  be  made 
available  to  the  image  coder.  It  is  the  finding  of  this  displacement 
vector  that  is  new  and  original  in  this  work.  Previous  work  only 
required  the  displacement  to  be  found  to  the  nearest  integer  multiple  of 
the  pixel  spacing.  This  may  suffice  for  some  applications,  but  where 
the  human  observer  is  the  final  link  in  the  coding  system,  these  integer 
only  displacements  tend  to  be  somewhat  annoying.  For  this  reason,  and 
others  concerned  with  actual  displacement  measurement,  the  noninteger 
portion  of  the  displacement  is  needed.  The  methods  of  pel  recursive  and 
coefficient  recursive  displacement  estimation  use  an  iterative  recursive 
minimization  procedure  to  converge  to  a  displacement  that  is  of 
noninteger  order. 

In  this  chapter  a  new  method  for  finding  the  displacement  vector  for 
use  in  a  motion  compensated  image  coding  system  is  introduced  and 
investigated.  The  method  of  prediction  coefficient  energy  concentration 
differs  from  the  earlier  recursive  methods  in  that  a  solution  for  the 

displacement  is  explicitly  defined  and  does  not  have  to  iterate  to  a 

solution.  Further,  it  differs  by  the  fact  that  no  actual  noninteger 
portion  of  the  displacement  actually  has  to  be  found.  Instead  this  is 
replaced  by  a  coefficient  prediction  process  that  obtains  the  same 

results,  or  in  many  cases  better,  with  no  increase  in  the  algorithm 

complexity. 
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metric  representing  each  corner  of  each  quadrant  where  the  center  value, 
being  common  to  all  quadrants,  may  be  neglected  as  is  shown  in  figure 
3.5.  Each  quadrant  then  is  represented  by  a  sum  of  three  metric  values 
and  the  minimum  quadrant  is  defined  to  be  the  one  with  the  minimum  sum. 
Ihis  in  essence  determines  the  integer  bounds  or  integer  displacement 
cell  for  the  noninteger  portion  of  the  displacement.  In  other  words, 
the  true  value  of  the  displacement  is  not  known,  but  it  is  assumed  to 
fall  within  the  bounds  of  the  displacement  cell.  Hence  the  locations  of 
the  four  delta  functions  are  defined  to  be  at  each  corner  of  the 
displacement  cell.  This  yields  the  locations  of  the  delta  functions, 
but  not  their  magnitudes.  The  method  for  determining  their  magnitudes 
is  given  in  the  following  section. 

At  this  point,  it  is  informative  to  look  at  a  two-dimensional  example 
for  the  displacement  of  a  single  pixel.  Figure  3.6(a)  presents  a 
reference  pixel  from  the  current  block.  Figure  3.6(b)  shows  the 
displacement  associated  with  this  pixel  as  well  as  the  integer  estimates 
for  the  X  and  T  displacements.  The  four  pixels  that  define  the  box  are 
used  in  a  linear  combination  to  estimate  the  reference  pixel. 
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FIGURE  3.6(a)  Current  Frame  Pixel  Location 
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FIGURE  3.6(b)  Displacement  Measurements  for  Reference  Pixel 

3  .5  PBKnrr-mg  COEFFICIENT  flUMBBATTON 

Having  bounded  the  displacement  estimation  in  both  directions,  the 
time-modified,  translation-updated  regression  model  is  given  as  before. 

1  1 

I(i,j,t)  -  £  £  a(m,n)I(i-K+m,  j-Lfn,  t—t)  +  e(i,j)  (3.35) 

a^O  a-0 

Where,  as  before,  K  and  L  designate  the  lover  bounds  on  the  displacement 
estimation.  As  vss  noted  earlier,  the  assumption  is  made  that  the 
displacement  remains  constant  over  the  block  and  hence  the  regression 
coefficients  vill  also  remain  constant  over  that  block. 

3.5.1  jtta  touiiiiaa 

Hie  solution  for  a  single  dimensional  regression  process  is  straight 
forward  and  vill  not  be  addressed  here.  Instead,  only  its  relationship 
to  the  problem  at  hand  vill  be  discussed. 

The  matrix  equation  normally  used  for  single  dimensional  regression 
can  be  written  as  follows 

Y  -  ZB  +  e  (3.36) 

vhere  Y  is  a  dependent  variable  vector  of  size  n.  X  is  the  sugmented 
independent  variable  matrix  of  size  n  by  p,  vhere  p  is  the  number  of 
regression  coefficients.  B  is  then  the  parameter  vector  or  vector 
containing  the  predictor  coefficients  and  is  also  of  size  p.  Finally  e 
is  the  error  vector  or  residual.  The  matrix  normal  equation  is  given  by 
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xxxb  »  xAr. 


(3 .37) 


The  least  squares  estimate  for  B  is  b  and  can  be  obtained  from 


(iTxr1xTr. 


(3.38) 


A  farther  restriction  is  placed  on  the  cnrrent  regression  problem  in 
that  it  is  an  antoregressive  series,  or  actually  an  autoregressive 
spatial-time  series,  and  current  values  are  predicted  from  past  time 
values  from  the  same  general  spatial  location.  The  autoregressive 
equation  can  be  written  as 


T  »  ZB  +  e. 


(3.39) 


Z  represents  a  shifted  version  of  7,  in  this  case  both  spatially  and 
temporally.  In  the  particular  case  at  hand,  the  derivation  has  to  be 
carried  one  step  further,  because  it  must  take  into  account  that  the 
data  is  two-dimensional  and  the  prediction  is  made  over  a  time  boundary. 
Note  that  the  regression  equation  given  in  (3.35)  is  not  in  matrix  form. 
Some  data  manipulation  is  required  if  standard  methods  for  regression 
analysis  are  to  be  used. 

3.5.2  Method  for  Data  Manipulation 

First  the  coefficient  matrix  a(i,j)  needs  to  be  placed  in  a  vector 
fitting  the  description  of  the  B  vector  in  (3.36),  therefore 


a(0.0) 

B  *  a(0.1)  (3.40) 

a(1.0) 

»(U)_ 

where  it  is  augmented  by  a(0) ,  the  intercept  or  bias  term.  Next,  the  Z 
matrix  needs  to  be  set  up  to  perform  the  shifting  that  is  accomplished 
by  the  double  summation  in  (3.35).  In  the  equation  for  the  Z  matrix 
that  follows,  the  t-x  factor  and  lower  integer  bound  terms  are  neglected 
for  simplicity. 
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I(i.j)  I<i.j+1)  I ( i+1 ,  j )  I( i+1, j+1) 

Ki.j+l)  I(i, j+2)  I(i+l,j+l)  I( i+1, j+2) 


I(i.j+N)  I(i,j+N+1)  I(  i+1 , j+N)  I(i+l.j+NH)  (3.41) 

I( i+1, j)  1(1+1, j+1)  I ( i+2 , j )  I(i+2, j+1) 


[1  I(i+N.j+N)  I( i+N,  j+N+1)  I(i+N+1,  j+N)  I(i+N+1.  j+N+1)  J 

Figure  3.7  is  provided  to  show  how  the  Z  matrix  is  generated  from  the 
two  dimensional  image  intensity  values.  Note,  the  first  column  of  Z  is 
augmented  with  l's  to  correspond  to  the  intercept  or  bias  term  of  the  B 
vector. 
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of  the  scan  diagrams  given  in  figure  3.7.  Figure  3.7(a)  represents  the 
pixel  groupings  by  row  for  the  intensities  from  the  previous  frame. 
Each  row  of  Z  is  represented  by  a  box  in  figure  3.7(a).  An  expanded 
version  of  one  of  these  boxes  is  shown  in  figure  3.7(b)  where  each 
corner  of  the  box  designates  a  column  in  the  Z  matrix.  Hence,  for  each 
row  in  Z  there  is  a  corresponding  group  of  four  pixels  whose  intensity 
values  are  to  be  placed  into  the  column  associated  with  the  number  given 
for  each  corner,  as  is  shown  in  figure  3.7(b).  Note  that  each  pixel 
intensity  may  be  used  up  to  four  times  in  the  Z  matrix.  Finally,  row 
scan  the  current  frame  block  placing  each  value  into  the  column  vector  T 
as  shown. 


I(i, j, t) 
1(1, j+l,t) 


I(i.N. t)  (3.42) 

I( i+1, j , t) 


[l(i+N,j+N.t)J 

The  scanning  diagram  for  the  current  frame  block  is  given  below  in 
figure  3.8.  The  residual  vector  e  is  defined  identically  to  that  of  the 
T  vector. 

With  each  of  the  variables  T,  Z,  B,  and  e  defined,  the  current 
problem  reverts  to  that  of  a  normal  one  dimensional  autoregression 
problem  involving  five  coefficients  and  hence  can  be  solved  as  such. 
The  least  squares  estimate  for  B  is  b  and  is  given  by 


(ZTZ)-1ZTY. 


(3.43) 


As  with  any  system  that  requires  a  matrix  inverse,  it  is  possible 
that  the  system  may  become  ill  conditioned  and  the  inverse  may  not 
exist.  In  the  current  system,  this  is  remedied  by  using  only  the 
estimate  for  the  integer  portion  of  the  displacement.  When  this  occurs, 
the  b  vector  is  set  to  an  identity  transfer  function,  that  is, 


FIGURE  3.8  T  Vector  Scan  Diagram 

bT  -  [0.1. 0.0,0]  (3.44) 

and  amounts  to  nalng  only  the  information  from  the  integer  displacement 
estimate  in  the  prediction  process. 

3 .6  BLOC!  mnafTPrraH  AND  BBSTniur  GENERATION 

After  the  completion  of  the  identification  portion  of  the  system,  it 
is  necessary  to  generate  the  data  required  for  transmission.  The  data 
required  by  the  receiver  is  broken  up  into  two  separate  parts.  The 
first  part  is  the  receiver  control  block,  which  contains  the  information 
generated  by  the  identification  portion  of  the  system.  The  remaining 
part  is  what  is  termed  the  residual  data  block.  It  contains  the 
information  required  by  the  receiver  to  update  or  correct  the  block 
estimate  when  based  only  upon  the  control  block  information.  The 
receiver,  as  well  as  the  receiver  portion  of  the  transmitter,  takes  the 
control  block  information  and  produces  an  estimate  of  the  current  block 
using  a  similarly  spatially  located  block  from  the  previous 
reconstructed  frame. 

1  1 

e(i.j.t)  *  I(i,j,t)  -  V  T  a(m,n)I( i-K+m, j-L+n, t~c)  -  a(0) 

opTo  n«0  (3.45) 

This  residual  term  is  then  quantized  for  transmission. 

In  order  to  fully  exploit  the  between  frame  redundancy,  checks  are 
placed  in  the  system  to  flag  the  types  of  changes  that  occur.  The  first 


check  determines  if  eny  motion  has  occurred  in  the  current  block.  If  no 
motion  is  found,  the  process  simply  goes  to  the  next  block  with  no 
required  data  transmission.  If  motion  has  occurred,  it  is  tested  to 
determine  if  an  integer  estimate  is  a  sufficiently  good  estimate.  If  it 
is  deemed  an  integer  only  displacement,  the  block  address  and 
displacement  are  transmitted  to  the  receiver.  When  integer  displacement 
is  not  accurate  enough,  the  prediction  coefficients  must  be  calculated. 
At  this  point  the  transmitter  must  send  the  block  address,  integer 
displacement,  quadrant,  predictor  coefficients  and  possibly  a  quantized 
version  of  the  prediction  error.  A  final  check  is  performed  to 
determine  if  the  error  is  still  too  large.  If  it  is,  the  primary 
blocksize  is  cut  in  half  and  the  process  is  retried  for  each  of  the  four 
sub-blocks. 

3.7  CODING 

As  noted  in  the  previous  section,  there  are  two  types  of  information 
that  must  be  transmitted  to  the  receiver.  The  first  is  the  control 
block  information,  which  contains  the  information  required  to  make  the 
motion  compensated  estimate  and  the  second  is  the  residual  data  block. 
The  residual  data  block  is  the  quantized  error  of  the  actual  prediction, 
along  with  the  coding  information  needed  by  the  quantizer.  Each  of 
these  will  now  be  discussed  in  detail  below. 

3.7.1  Control  Blook  Information 

The  control  block  information  data  sequence  contains  the  set  up 
variables  needed  by  the  receiver  in  order  to  start  the  prediction 
process.  The  data  sequence  is  of  variable  length,  depending  on  the  mode 
of  operation,  of  which  there  are  four.  As  noted  earlier,  if  no  change 
has  occurred  for  a  block,  then  the  receiver  requires  no  update. 

The  four  different  modes  of  operation  can  be  separated  into  two 
groups  of  two  categories.  The  two  groups  are  integer  and  noninteger 
displacements,  while  the  two  categories  are  the  different  blocksizes  of 
8  by  8  and  16  by  16.  Table  3.1  shows  the  bit  requirements  for  each  mode 
of  operation  and  the  bit  breakdown.  The  assignment  of  the  first  10  bits 
remain  the  same  for  all  four  modes.  These  first  10  bits  consist  of  a 
bit  for  integer/noninteger  mode,  a  blocksize  bit  for  blocksize 


determination  and  4  bits  each  for  the  integer  X  and  7  displacement. 
This  will  allow  for  a  shift  in  either  direction  of  from  minus  seven  to 
pins  seven  pixels.  The  remainder  of  the  bit  structure  differs  in 
accordance  with  the  blocksize  bit  and  integer  displacement  bit.  These 
remaining  bits  contain  the  block  address  for  the  current  block  and,  if 
the  noninteger  bit  is  set,  the  quadrant  number.  Four  of  the  predictor 
coefficients  are  coded  into  this  section,  while  the  fifth,  since  it  is 
used  by  the  residual  quantizer,  is  coded  in  the  residual  block  section. 
The  number  of  bits  used  is  based  on  an  image  size  of  256  pixels  square 
or  smaller.  If  a  larger  image  is  used,  then  appropriate  changes  will 
require  the  enlargement  of  the  block  address. 

This  control  block  is  essentially  the  overhead  required  for  the 
prediction  process  and  in  itself  produces  a  very  low  data  rate.  This 
overhead  is  zero  for  the  motionless  blocks  and  can  range  from  about  0.07 
up  to  0.98  bits  per  pixel  for  blocks  that  contain  motion.  Hence,  the 
predictor  overhead  will  normally  contribute  only  a  small  portion  of  the 
total  data  rate  requirements.  The  majority  of  the  bandwidth  required  is 
a  result  of  the  prediction  error  and  hence  is  required  by  the  residual 
data  block. 


TABLE  1 


INTEGER 


NON-INTEGER 


3-Integer  only  0-Integer  only  1-Non-Integer  1-Non-Integer 
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3.7.2  Residual  Data  Block 


Then  the  euaulative  error  between  the  predicted  block  end  actual 
block  exceeds  a  set  threshold,  a  residual  block  Bust  also  be 
transaitted.  This  data  is  used  to  correct  some  of  the  errors  that  occur 
in  the  prediction  process.  The  block  consists  of  two  parts.  The  first 
contains  the  Bean  and  variance  of  the  residual  aignal  and  the  nuaber  of 
quantizer  levels.  The  reaaining  data  is  the  quantized  residual  data. 
Fixed  rate  quantizers  are  inappropriate  because  the  residual  signal  is 
noir-stationary .  A  variable  rate  quantizer,  based  on  the  aean  and 

variance  of  the  residual  signal  and  the  predictor  gain,  perforaed  quite 
well  in  spite  of  the  non-stationary  signal.  For  very  saall  variance 
errors  no  residual  was  needed.  For  larger  errors,  up  to  7  bits  were 
available  per  pixel. 

Due  to  the  very  large  variation  in  the  quality  of  the  predicted 

signal,  an  adaptivr-variable-length  quantizer  is  required  in  order  to 
aaintain  a  ainiaal  data  rate  for  a  specified  overall  picture  quality. 
Here,  as  in  auch  of  the  other  work  in  iaage  processing,  there  is  no  good 
clear  cut  uuaerical  indication  of  picture  quality.  For  the  work 

reported  here  two  statistical  indicators  were  jointly  used  for  internal 
judgeaent  of  picture  quality.  It  is  upon  the  basis  of  this  judgeaent 

that  the  residual  data  rate  is  deterained.  The  first  statistical 

indicator  is  the  error  signal  variance.  It  is  expected  that  a  low 
variance  value  aeans  good  iaage  reproduction  whi!e  a  large  error 
variance  indicates  a  large  prediction  error  and  hence  a  large  residual 
data  rate  requireaent.  The  second  statistical  indicator  is  the 

prediction  gain  or  predictor  signal  to  noise  ratio.  Thile  this  does 

take  into  account  the  error  variance  froa  above,  it  also  takes  into 
account  the  signal  variance  and  is  defined  by 

GAIN  -  10Log(var(signal) /var(error) ) .  (3.46) 

Using  these  two  quality  indicators  as  fidelity  criteria,  the  adaptive 
quantizer  deteraines  the  required  nuaber  and  placeaent  of  quantizer 
output  levels  in  order  to  achieve  a  specified  output  signal  to  noise 
ratio.  The  residual  variance  is  used  as  a  simple  yes/no  indicator  for 
the  quantizer.  That  is,  if  the  error  variance  exceeds  the  preset 
threshold,  the  quantizer  will  be  used  to  transait  the  residual.  The 
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predictor  signal  to  noise  ratio  is  then  used  by  the  quantizer,  along 
with  the  minimum  acceptable  signal  to  noise  ratio,  to  determine  how  much 
the  predicted  signal  needs  to  be  improved.  This  is  equivalent  to 
determining  the  quantizer  level  requirements.  The  number  of  levels  is 
determined  by  finding  the  minimum  number  of  levels  needed  to  obtain  the 
preset  output  signal  to  noise  ratio,  provided  the  error  distribution  is 
Gaussian.  Even  though  the  global  error  distribution  may  be  gaussian, 
rarely  will  the  local  distribution  be  so.  For  this  reason,  rarely  will 
the  output  signal  to  noise  ratio  match  that  which  was  preset.  However, 
this  will  bound  the  actual  data  rate  requirements. 

The  quantizer  employed  in  the  actual  implementation  is  a  multiple 
level  and  hence  variable  rate.  The  quantizer  chooses  the  thresholds  and 
output  levels  optimally  based  on  a  distortion  measure  defined  to  be  the 
sum  of  the  absolute  values  of  the  error.  The  choice  of  this  distortion 
measure  over  the  standard  mean  squared  error  results  in  a  decrease  of 
the  granular  noise  in  the  areas  of  the  image  with  small  signal  variance. 
Figures  3.9(a)  and  (b)  are  provided  to  show  the  differences  in  output 
error  distributions  for  the  different  distortion  measures.  Each  is 
assumed  to  have  a  Normal (0,1)  input  error  distribution.  The  curves 
labelled  0,1, 2, 3,  and  4  bits  are  the  error  output  distributions  after 
modification  by  the  selected  quantizer.  If  an  infinite  number  of  bits 
were  used,  the  output  error  distribution  would  be  a  delta  function 
centered  about  0.  Note  that  the  mean  square  method  tends  to  minimize 
the  area  under  the  tails  of  the  distribution,  while  the  absolute  value 
method  tends  to  concentrate  more  on  the  center  portion  of  the 
distribution.  As  noted  earlier,  the  number  of  quantizer  levels  required 
is  based  upon  the  error  variance  and  the  predictor  signal  to  noise 
ratio. 

When  the  error  variance  exceeds  the  set  threshold,  the  difference 
between  the  requested  gain  and  the  actual  predictor  gain  is  in  a  sense  a 
'gain'  that  must  be  generated  by  the  quantizer.  Gain  is  a  term  not 
normally  associated  with  quantizers,  but  in  this  case  it  is  used  to 
compare  the  output  quality  with  a  given  quantizer  to  the  quality  that 
would  be  present  if  no  quantizer  were  used.  That  is  to  say,  how  much  is 
the  signal  improved  by  using  a  quantizer  to  send  the  error  ss  compared 
to  not  sending  the  error  at  all. 


This  gain  is  determined  ss  s  function  of  the  number  of  output  levels. 
It  is  then  s  simple  procedure  to  determine  how  many  levels  are  required 
to  meet  gain  requirements  to  obtain  this  reproduced  image  fidelity 
criterion.  A  plot  of  this  gain  as  a  function  of  output  levels  is  given 
in  figure  3.10.  From  this,  for  a  required  gain,  the  number  of  output 
levels  can  be  determined.  The  similarity  of  this  curve  to  the  log  rate 
distortion  curve  is  apparent. 

3 


FIGURE  3.10  Quantizer  Gain  as  a  Function  of  Number  of  Output  Levels 
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Chapter  IV 

EXPERIMENTAL  OPERATION 


The  presentation  of  results  obtained  from  most  forms  of  image 
processing  is  a  difficult  one,  but  even  more  so  when  the  data  consists 
of  tine  varying  image  sequences,  as  is  the  case  in  motion  compensated 
image  coding.  Comparisons  among  published  results  are  even  more 
difficult  because  of  the  lack  of  standard  presentation  procedures  and 
the  lack  of  standard  test  image  sequences.  The  performance  of  any 
motion  compensated  coding  system  is  extremely  scene  or  sequence 
dependent  and  hence,  this  dependence  as  a  function  of  system  and  scene 
characteristics  will  be  presented. 

The  image  sequences  used  for  this  system  performance  analysis  consist 
of  41  images  of  size  128  by  192  pixels  quantized  to  8  bits.  The  first 
sequence  is  termed  'SLOW  PHONE'  and  consists  of  a  woman  talking  on  a 
phone  while  slowly  rotating  and  moving  the  phone.  The  second  sequence, 
termed  'FAST  PHONE',  is  similar  to  the  first  but  contains  faster  and 
more  motion.  Both  sequences  contain  very  complex  motion  with 
foreground-background  interaction  and  are  similar  to  those  that  would  be 
encountered  in  a  video-phone  setting. 

A  block  diagram  showing  each  of  the  subsystems  mentioned  in  the 
previous  chapter  and  their  relationship  to  one  another  is  given  in 
figure  4.1.  Note  that  the  block  diagram  presents  each  subsystem 
serially  connected  with  no  indication  for  mode  of  operation.  The  mode, 
recall  there  are  four,  is  controlled  by  the  data  control  and  logic 
block.  Note  also  that  the  system  can  be  easily  implemented  in  parallel 
in  terms  of  block  operations  because  one  block  prediction  is  not  based 
on  the  previous  block  from  the  same  frame,  but  only  on  blocks  from  the 
previous  reconstructed  frame. 

As  stated  above,  the  presentation  of  the  results  for  time  varying 
imagery  work  is  very  difficult.  When  images  from  the  sequence  are 
viewed  singularly,  they  appear  'cleaner'  or  less  noisy  than  they  would 
if  viewed  in  the  actual  sequence.  Also,  some  prediction  errors  that  do 
appear  in  the  single  frames  are  not  as  noticeable  when  the  frames  are 
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DIGITAL  CHANNEL 


BLOCK 


Upon  return  from  the  metric  subroutine,  sn  estimate  of  the  integer 
displacement  is  provided.  The  subroutine  also  supplies  a  variable  that 
is  equal  to  the  value  of  the  metric  at  its  minimum  location.  This  value 
can  be  used  to  determine  how  good  of  an  estimate  of  the  current  block 
can  be  generated  using  the  previous  frame  and  the  integer  displacement 
estimate.  For  many  blocks  this  estimate  may  prove  to  be  sufficient. 
For  others,  it  indicates  that  more  information  will  be  required  in  the 
image  reconstruction.  For  the  blocks  that  the  estimate  is  sufficient, 
no  further  displacement  estimation  is  required.  For  the  others,  further 
coeiputations  are  required  to  improve  the  current  frsme  estimate.  This 
further  computation  may  involve  determining  a  better  estimate  for  the 
displacement  vector  or,  as  in  the  method  presented,  solve  for  the 
regression  parameters. 

When  an  integer-only  displacement  estimate  is  not  good  enough,  the 
next  step  is  to  calculate  the  best  estimate  for  the  quadrant  in  which 
the  true  displacement  falls.  This  is  accomplished  through  the  use  of 
the  LOCATE  subroutine  as  given  in  appendix  B.  The  accompanying 
flowchart  is  provided  in  figure  4.7.  Recall  from  chapter  3  that  this 
quadrant  is  generated  by  sususing  the  values  of  the  metric  that  surround 
the  best  integer  estimate.  The  quadrant  is  used  to  bound  the 
displacement  estimate  to  a  single  integer  displacement  cell.  This 
subroutine  locates  the  quadrant  with  the  lowest  sum  of  surrounding  pixel 
values.  It  does  not  determine  the  values  for  the  coefficients. 
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FIGURE  4.6(d)  Plot  of  Actntl  Metric 


estimate  of  the  displacement.  This  presentation  is  achieved  with  the 
ose  of  subroutine  INVERSE  and  MPLOT  as  given  in  appendix  B.  Figures  4.6 
(a)  through  (d)  are  examples  of  various  actual  metrics.  From  these,  the 
geometrical  nature  of  the  metric  function  can  be  seen. 

The  plots  are  arranged  such  that  the  base  of  each  plot  represents  the 
value  of  the  integer  shift.  Each  base  axis  represents  a  direction  of 
shift.  The  values  of  shift  for  these  examples  range  from  minus  nine  to 
plus  nine  in  both  the  X  and  Y  direction.  The  magnitude  of  the  plot  can 
be  thought  of  as  a  measure  of  similarity  between  the  current  block  and 
various  shifted  blocks  from  the  previous  frame. 


difficult  in  practice.  Following  the  flow  chart  of  figure  4.5,  the 
process  starts  with  setting  an  initial  guess  for  the  displacement  vector 
estimate.  This  is  not  used  as  starting  point  for  iteration,  but  rather 
as  a  method  to  save  time.  This  value  can  be  set  to  zero,  if  no  motion 
is  expected,  or  to  the  value  of  the  previous  block  estimate.  The 

subroutine  will  test  the  hypothesis  that  this  is  the  correct  integer 
displacement,  and  if  it  falls  within  the  threshold  bounds,  will  exit  the 
subroutine  setting  the  actual  estimate  to  the  initial  guess.  If  the 
metric  value  falls  outside  the  threshold,  the  entire  procedure  must  be 
executed  for  all  allowable  possibilities  of  block  shift.  The  procedure 
is  executed  as  follows.  First  the  shift  vector  estimate  must  be  set  to 
the  lowest  possible  value  in  both  the  X  and  Y  directions.  Solve  for  the 
metric  or  risk  function  based  on  this  shift  vector.  Recall  from  Chapter 
3  that  the  metric  is  defined  to  be  a  function  of  the  sum  of  the 

thresholded  differences  between  the  current  frame  block  and  the  previous 
frame  pixel  neighborhood.  The  shift  vector  is  then  incremented  and  a 
new  metric  value  calculated.  The  procedure  is  continued  until  all 
allowable  pixel  displacements  have  been  tried.  At  this  point,  a  matrix 
of  metric  values.  whose  entry  position  defines  the  associated 

displacement  estimate,  is  examined  for  a  minimum  value.  The  location  of 
the  minimum  value  designates  the  integer  displacement  estimate.  If 

there  is  more  than  one  minimum  value,  the  entire  metric  matrix  is 
recalculated  using  a  smaller  allowable  error  threshold.  If  after  this 
is  tried,  and  the  allowable  error  cannot  be  decreased  and  there  are 
multiple  minimum  points,  the  smallest  displacement  value  is  chosen  to  be 
the  actual  estimate.  Along  with  the  location  of  the  smallest  metric 
value,  the  actual  value  of  the  metric  is  also  returned  to  the  calling 
program  as  an  aid  in  quality  measurement.  The  program  is  now  ready  to 
return  to  the  main  driver  program  and  continue. 

Before  continuing  with  the  program  procedure,  it  is  informative  to 
look  at  some  examples  of  the  metric  matrices  generated  by  the  algorithm. 
Recall  that  the  actual  implementation  looks  for  a  minimum  value  of  the 
metric  matrix,  but  for  ease  of  presentation  the  metric  matrix  has  been 
modified.  The  modification  is  accomplished  by  reversing  the  magnitude 
order  of  the  data  and  rescaling  to  fit  the  0  to  1  range.  The  reverse  in 
order  exchanges  places  of  the  minima  and  maxima.  Hence,  in  the 
presentation  the  maximum  value  of  the  function  defines  the  best  integer 


FIGURE  4.5  Metric  Subroutine  Flow  Chart 
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Following  constant  initialization,  the  single  image  is  broken  up  into 
many  small  sub-images  of  size  16  by  16  and  blockwise  scanned  top-to- 
bottom  and  lef t-to-right.  The  following  discussion  refers  to  operations 
on  a  single  block. 

First  it  must  be  determined  if  the  particular  block  is  a  border 
block.  That  is,  does  it  border  the  outside  of  the  image?  The  reason 
for  this  is  that  border  blocks  request  data  from  the  previous 
reconstructed  pixel  neighborhood  that  do  not  exist.  This  is  dictated  by 
the  possibility  of  a  plus  or  minus  p  pixel  shift  in  both  directions. 
Bather  than  setting  the  nonavailable  values  to  zero,  it  has  proved 
advantageous  to  set  these  terms  to  the  mean  value  of  the  remaining  pixel 
neighborhood. 

4.3  METRIC  BBMRBATTOM 

Given  the  current  block  and  the  previous  pixel  neighborhood,  the  next 
step  is  to  calculate  the  difference  metric.  The  flow  chart  for 
subroutine  METRIC  is  provided  in  figure  4.3  while  the  subroutine  listing 
is  provided  in  appendix  B. 

Subroutine  METRIC  serves  a  dual  purpose.  It  is  used  to  determine  the 
best  estimate  for  the  integer  portion  of  the  displacement  based  upon  the 
arithmetic  metric  used.  It  is  also  used  to  determine  if  this  integer 
estimate  is  a  sufficiently  good  estimate.  In  the  discussion  to  follow 
the  absolute  value  threshold  method,  or  MAP  estimator  given  in  section 
3.3.3,  will  be  used.  In  the  subroutine,  like  the  main  driver,  it  is 
important  to  treat  border  and  nonborder  blocks  differently.  The 
execution  is  similar  for  both  border  and  nonborder  blocks,  with  the 
exception  of  a  correction  factor  that  takes  into  account  the  possibility 
of  a  smaller  number  of  terms  used  in  the  border  block  calculations.  The 
correction  factor  is  used  to  weight  the  estimates  based  on  fewer  terms 
differently  than  it  would  if  the  entire  previous  data  field  had  been 
available. 

It  is  this  portion  of  the  algorithm  that  is  the  most  calculation 
intensive,  and  hence  the  slowest  to  execute.  Here  again  the  algorithm, 
although  implemented  in  series,  lends  itself  to  parallel  operations. 
The  implementation  of  the  algorithm  is  simple  in  concept,  but  more 


FIGURE  4.4  Main  Control  Program  Flow  Chart 


This  method  for  motion  compensated  image  coding  lends  itself  very 
well  to  modular  programming  design.  The  system  must  perform  serially  in 
terms  of  most  within  block  operations,  but  is  easily  adaptable  to 
perform  parallel  or  concurrent  operations  across  blocks.  Because  the 
algorithm  was  modelled  and  simulated  in  software,  only  serial  operations 
can  be  used  and  hence,  the  discussion  to  follow  will  be  restricted  to  a 
purely  serial  implementation.  Notes  will  indicate  the  parts  of  the 
implementation  that  are  easily  converted  to  parallel  operations. 

As  noted  earlier,  it  is  assumed  that  both  the  transmitter  and 
receiver  have  complete  knowledge  of  the  initial  image  and  motion 
compensated  image  coding  starts  with  the  second  frame  of  the  sequence. 
The  actual  FORTRAN  code  used  for  the  simulation  is  provided  in  appendix 
B  and  will  be  referred  to  throughout  this  chapter.  Flowcharts  for  both 
the  main  program  and  some  of  the  subroutines  will  be  provided  in  the 
text  as  required. 

4 .2  MAIN  CONTROL  PROGRAM 

Figure  4.4  provides  a  flow  chart  of  the  main  control  program  and 
relates  to  the  program  listing  MAIN  in  appendix  B.  The  program 
initialization  section  sets  up  the  required  run  time  constants  used  for 
both  variable  assignment  and  program  control.  Some  of  the  more 
important  constants  will  be  discussed  below. 

NPBITS  is  an  integer  number  used  to  designate  the  number  of  bits  used 
to  quantize  the  predictor  coefficients.  The  quantization  method  assumes 
a  uniform  distribution  from  -2  to  2.  The  actual  number  of  quantizer 
bits  is  one  greater  than  NH3ITS  to  allow  for  a  sign  bit.  The  variables 
VAREST,  STDERR,  and  AVGERR  are  set  constants  used  as  threshold  values 
throughout  the  program.  VAREST  is  an  estimate  of  the  average  difference 
between  two  identical  frames  when  viewed  at  the  output  of  the  imaging 
device,  or  in  other  words  an  estimate  of  the  imager  noise.  STDERR  and 
AVGERR  are  threshold  values  for  the  error  terms.  ISIZE  and  JSIZE  are 
simply  used  to  identify  the  image  size  and  may  be  changed  to  fit  any 
image  size  as  long  as  the  DIMENSION  statement  settings  are  likewise 
adj  usted. 
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viewed  in  sequence  because  of  the  integrating  ability  of  tbe  eye.  For 
simple  comparison  purposes,  six  frames  from  both  of  tbe  original 
sequences  will  be  presented.  Figures  4.2(a)  through  (f)  are  frames  from 
the  'SLOW  PHONE'  sequence  and  4.3(a)  through  (f)  are  from  the  'FAST 
PHONE'  sequence.  It  is  important  to  note  that  much  of  the  image  quality 
is  lost  in  the  photographic  process  and  the  associated  copying  and 
printing. 

The  theoretical  system  presented  in  Chapter  3  has  been  programmed  and 
simulated  using  FORTRAN  and  hence  was  not  executed  in  hardware  or  real 
time.  More  on  the  actual  implementation  of  the  algorithm  will  be 
presented  in  the  later  sections. 

The  algorithm  simulation  starts  with  the  assumption  that  both  the 
transmitter  and  receiver  have  full  knowledge  of  the  first  frame.  The 
following  frame  starts  the  process  of  motion  compensated  image  coding. 
Transmission  of  the  first  frame  may  be  accomplished  using  normal  coding 
techniques  or  transmitted  using  ordinary  PCM  methods.  The  process  from 
this  point  onward  assumes  that  one  frame  is  read  at  a  time  and  processed 
before  the  next  frame  can  be  read.  Although  none  is  used  in  this 
simulation,  it  is  understood  that  in  a  real  world  implementation  a  fixed 
length  buffer  would  have  to  be  incorporated  into  the  hardware.  It  is 
assumed  that  an  infinite  length  buffer  is  available  in  the  simulation. 
Although  not  used  for  buffer  control,  a  very  good  built  in  feature  of 
the  implementation  is  the  setting  of  a  required  output  signal  to  noise 
ratio.  If  a  fixed  length  buffer  were  used,  the  varying  of  this  required 
SNR  could  function  as  a  buffer  control  parameter.  That  is,  when  the 
buffer  approaches  full  the  required  SNR  could  be  allowed  to  drop.  When 
the  buffer  neared  empty,  the  output  SNR  could  then  be  increased.  This 
would  nearly  eliminate  the  possibility  of  a  buffer  overflow,  and  hence  a 
loss  of  data. 


The  generation  of  the  delta  function  magnitudes  is  based  on  linear 
regression  theory  and  is  accomplished  in  the  subroutine  COEPGN.  The 
magnitudes  of  the  delta  functions  are  used  as  weighting  values  to 
estimate  the  current  frame  from  the  past  frame.  Given  the  integer 
displacement  value  and  the  magnitudes  of  each  of  the  delta  functions, 
the  system  is  able  to  make  a  good  prediction  of  the  current  block. 

Because  the  displacement  has  been  bounded,  the  four  pixels  that  are 
used  to  estimate  the  current  pixel  are  known.  What  is  not  known,  is  how 
these  pixel  intensities  are  combined  to  produce  the  current  estimate. 
The  weighting  values  define  how  the  previous  pixels  are  to  be  combined 
and  are  determined  by  solving  for  the  parameters  of  the  autoregressive 
spatial-time  series. 

4.6  riTBBKNT  BLOCI  ORWRHATIQN 

Subroutine  PRED  is  the  portion  of  the  program  that  is  used  to 
calculate  the  current  block  estimate.  Along  with  this  estimate,  an 
error  block  is  also  calculated.  From  the  error  block,  the  error  mean 

and  variance  must  be  determined.  The  mean  and  variance  are  used  as  an 

indication  of  how  well  the  prediction  process  performed.  If  both  the 
mean  and  variance  are  within  set  limits,  the  process  will  aimply  go  on 
to  the  next  block.  If  the  error  is  deemed  too  large,  the  entire  process 
up  to  this  point  is  redone  on  a  smaller  blocksize  of  8  by  8 .  This  says 
that  the  single  Id  by  Id  block  is  broken  up  into  4  sub-blocks  and  the 
algorithm  is  executed  for  each  of  the  sub-blocks.  If  the  blocksize  is 
already  8  by  8  and  the  error  is  still  too  large,  the  process  must 
generate  a  quantized  error  or  quantized  residual  term. 

The  error  in  prediction  can  come  about  from  many  sources,  some  of 

which  are  listed  below.  The  motion  that  occurs  may  not  fall  into  the 
category  of  pure  translation.  That  is,  the  motion  may  consist  of 

various  types  of  rotation.  The  objects  could  have  rotated  in  the  plane 
of  view  or  in  any  of  the  planes  perpendicular  to  the  plane  of  view, 
hence  taking  into  account  the  3-D  aspect  of  the  objects.  The  model 
cannot  account  for  rotation.  Hence,  when  this  occurs  errors  will  result 
that  must  be  quantized  and  transmitted.  The  other  major  source  of  error 


generation  is  from  the  problem  of  background/ foreground  interaction. 
Recall  that  the  motion  model  was  developed  for  the  pure  translation,  no 
background  case  and  hence  only  optimally  estimates  this  type  of  motion. 
When  the  actual  data  does  not  fit  the  model,  it  may  be  necessary  to 
quantize  and  transmit  the  error. 


4.7  ERROR  QUANTIZATION 

The  error  term  resulting  from  the  subtraction  of  the  predicted  block 
from  the  actual  block  is  very  scene  dependent.  Because  of  this  the 
statistics  for  the  error  block  are  widely  varied.  The  quantizer  used 
for  the  residual  encoding  is  very  important  because  it  determines,  to  a 
large  extent,  the  overall  system  data  rate.  The  variance  of  the  error 
is  a  good  indicator  of  image  quality.  That  is,  if  the  variance  is  low 
the  image  reproduction  should  be  good.  On  the  other  hand,  if  the 
variance  of  the  error  is  large,  it  is  expected  that  the  reproduced  image 
will  be  somewhat  degraded.  It  only  seems  logical,  from  the  above 
discussion,  that  the  quantizer  used  should  be  able  to  adapt  to  the 
changes  in  statistics  of  the  error  and  transmit  only  the  data  reqnired 
to  reach  some  fidelity  criterion. 

Subroutines  QOANTI  and  Q0AN1Z  are  used  in  the  quantization  of  the 
error  signal.  QOANTI  determines  the  amount  of  gain  the  quantizer 
produces  for  the  various  number  of  output  levels.  QOANTZ  performs  the 
actual  quantization  and  coding  of  the  residual. 


Chapter  V 

RESULTS  AND  CONCLUSIONS 

The  presentation  of  results  from  time  sequential  image  processing, 
and  even  image  processing  in  general,  proves  to  be  a  difficult  task. 
When  a  human  viewer  is  the  final  link  in  this  imaging  system,  the 
important  criterion  is  not  a  set  of  numbers  calculated  from  various 
system  parameters,  but  rather  the  important  criteria  is  how  does  the 
viewer  subjectively  judge  the  quality  of  the  final  output.  Although 
there  appears  to  be  some  correlation  between  some  distortion  measures, 
such  as  the  mean  square  error,  and  the  subjective  quality,  this 
correlation  is  not  very  strong.  Another  way  to  put  this  is,  a  decrease 
in  the  mean  square  error  may  not  always  mean  an  increase  in  judged 
subjective  quality. 

A  presentation  of  output  results  of  single  images,  not  taken  from  an 
image  sequence,  poses  a  simpler  problem  than  those  taken  from  image 
sequences.  Each  image  can  be  viewed  subjectively  apart  from  all  others 
because  they  are  not  dependent  on  images  that  come  both  before  and  after 
in  time.  This  method  of  presentation  does  not  hold,  as  well  as  is  not 
feasible,  for  sequential  image  data.  When  the  sequences  are  composed  of 
thirty  or  more  frames  per  second,  the  amount  of  space  required  to 
present  them  in  itself  creates  a  problem.  But  more  importantly,  the  eye 
and  associated  biological  processors  prefer  that  the  images  be  presented 
in  a  sequential  manner  in  the  same  spatial  location.  It  appears  that 
the  integrating  ability  of  the  eye  plays  an  important  role  in  the 
subjective  image  judgement  of  time  sequential  imagery. 

The  output  images  will  be  presented  in  a  fashion  similar  to  the 
presentation  of  the  input  images  provided  in  chapter  four.  One  should 
note  that  quality  degradation  occurs  at  every  step  of  the  image  transfer 
process.  That  is,  exposure,  development,  printing  and  copying  all 
degrade  the  true  quality  of  the  actual  image. 

Even  though  distortion  measures  such  as  the  mean  square  error  may  not 
be  directly  related  to  the  subjective  quality  of  the  output  image,  it  is 
an  indicator  that  allows  for  a  comparison  among  different  runs  of  the 


same  program  with,  varied  parameters.  It  does  function  as  a  means  to 
compare  outputs  on  a  more  numerical  level.  It  seems  logical  that  if  the 
system  and  algorithm  remain  the  same,  with  only  minor  modifications  in 
some  system  parameters,  that  the  outputs  can  he  compared  on  a  more 
quantitative  basis.  Stating  this  simply,  the  image  sequence  with  the 
best  overall  output  signal  to  noise  ratio  should  be  the  image  sequence 
judged  best  in  a  subjective  analysis. 

Another  problem  in  the  presentation  of  results  arises  from  the  very 
large  variation  in  possible  input  image  sequences.  Because  the  data 
compression  that  occurs  in  the  system  is  a  direct  result  of  the 
exploitation  of  the  between  frame  correlation,  the  amount  of  data 
compression  that  can  occur  is  directly  related  to  how  well  adjacent 
frames  are  correlated.  It  has  been  shown  in  intraframe  image  coding 
techniques  [8],  that  the  majority  of  the  information  of  the  image  is 
contained  in  the  edges  of  the  image.  Along  with  the  majority  of  the 
information,  the  edges  also  require  the  majority  of  the  bandwidth  in  a 
data  compression  scheme.  The  same  can  be  said  for  interframe  image 
coding  if  one  defines  what  is  meant  by  'edge'  in  interframe  terms.  If 
an  intraframe  edge  is  defined  to  be  a  boundary  between  two  regions  of 
near  uniform  luminance  and  can  be  estimated  by  using  a  spatial 
derivative,  the  interframe  edge  can  be  defined  as  a  boundary  between 
spatial  regions  of  varying  intensity  temporally  separated  and  can  be 
estimated  with  a  temporal  derivative.  What  this  says  is  that  the  edge, 
in  interframe  terms,  occurs  at  the  portions  of  the  image  that  change 
between  frames.  This  idea  can  be  further  restricted  if  the  method  of 
interframe  image  coding  is  motion  compensated.  Now  the  edge  can  be 
defined  to  occur  at  the  regions  of  the  image  that  interface  the 
background  and  foreground.  Hence,  if  this  interframe  edge  were  viewed, 
the  edges  would  appear  at  points  where  the  moving  objects  and  background 
meet.  Just  as  in  the  intraframe  case,  the  temporal  edge  requires  the 
largest  amount  of  bandwidth  for  transmission,  and  hence  contains  the 
majority  of  the  new  information  contained  in  each  image. 

Because  the  amount  of  information  that  must  be  transmitted  is 
directly  related  to  the  quantity  of  temporal  edges  and  hence  motion,  the 
amount  of  motion  that  occurs  in  a  given  image  sequence  should  provide  a 
good  indication  of  the  data  rate  requirements  of  the  system. 


For  presentation  purposes,  only  a  fraction  of  the  total  number  of 
images  can  be  presented.  Figures  5.1(a)  through  (f)  are  a  selected  set 
of  output  images  from  the  'SLOW  PHONE'  sequence  with  the  following  set 
of  selected  parameters.  These  parameters  will  be  changed  to  try  and 
show  the  effects  of  the  various  thresholds  and  quality  settings  on  the 
quality  of  the  output. 


NFS ITS  -  10 
VAREST  -  1.0 
STDESR  “3.0 
AVGERR  -  3 .0 
SNESET  -  24.0 


Figures  5.2(a)  through  (f)  are  also  from  the  ' SLOW  PHONE'  sequence 
with  the  following  set  of  selected  parameters. 


NFBITS  -  10 
VAREST  -1.5 
STDESR  -4.0 
AVGERR  -4.0 
SNRSET  -  21.0 


Figures  5.3(a)  through  (f)  are  selected  images  from  the  'SLOT  PHONE' 
sequence  with  the  following  set  of  selected  parameters. 


NFBITS  -  10 
VAREST  -1.0 
STDERR  -  2.0 
AVGERR  -2.0 
SNRSET  -  27.0 


Figures  5.4(a)  through  (f)  are  also  images  from  the  'FAST  PHONE’ 
sequence  with  the  following  set  of  selected  parameters. 


NPBITS  -  10 
VAREST  -1.0 
STDERR  -  3 .0 
AVGERR  -  3 .0 
SNRSET  -  24. 


As  stated  earlier,  the  most  important  aspect  of  the  output  is  the 
subjective  quality  of  the  images,  but  other  parameters  of  the  system  may 
provide  insight  into  functioning  of  the  algorithm.  For  this  reason,  a 
set  of  eleven  plots  for  each  of  the  four  output  sequences  will  be 
presented.  The  plots  will  be  presented  in  groups  of  four,  each  one 
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plotting  the  same  parameter  aa  the  others  in  the  group,  but  with  a 
different  image  set  or  initial  paraaeter  setting. 

The  first  set  of  plots,  figures  5.5(a)  through  (d),  provide  a  measure 
of  the  amount  of  motion  that  is  present  in  the  image  sequence.  For  the 
purpose  of  the  plots,  motion  is  defined  to  have  occurred  within  a  block 
if  any  information  about  that  block  must  be  transmitted.  This  includes 
cases  of  integer-only  displacement,  nonrinteger  displacement  and 
possible  cases  where  no  motion  has  occurred,  but  a  residual  data 
sequence  is  required.  For  counting  purposes,  a  block  is  defined  as  an 
eight  by  eight  sub-image.  The  results  are  plotted  as  a  percentage  of 
the  maximum  number  of  sub-blocks  possible.  For  the  image  size  used, 
with  an  eight  by  eight  blocksize,  the  total  number  possible  is  384.  The 
average  number  of  blocks  considered  to  be  in  motion  in  the  sequence  of 
forty  images  is  printed  in  the  upper  right-hand  corner  of  each  figure. 
Here,  as  in  the  remaining  ten  plot  sets,  the  (a)  plot  relates  to  the 
output  figures  provided  in  figures  5.1(a)  through  (f) .  Plots  labelled 
(b),  (c),  and  (d)  then  relate  to  figures  5.2(a)  through  (f)  to  5.4(a) 
through  (f)  respectively. 

The  second  set  of  plots,  figures  5.6(a)  through  (d),  present  the 
instantaneous  data  rate  as  a  function  of  frame  number.  Like  the 
previous  set  of  plots  for  the  block  rate,  the  plots  start  from  a  small 
value  and  quickly  rise  to  a  more  stable  value.  This  can  be  partially 
explained  by  noting  that  it  may  take  multiple  frames  to  detect  and 
correct  for  very  small  amounts  of  motion  between  frames. 

The  next  set  of  plots,  figures  5.7(a)  through  (d),  simply  combine  the 
two  previous  sets  of  plots  in  such  a  way  so  that  the  interaction  of  the 
block  rate  and  data  rate  can  be  more  easily  seen.  As  would  be  expected, 
the  points  tend  to  cluster  about  a  line  angled  from  the  lower  left  to 
the  upper  right.  This  simply  shows  that  an  increase  in  motion  will 
require  an  increase  in  data  rate. 

Figures  5.8(a)  through  (d)  provide  a  numerical  indication  of 
instantaneous  system  distortion.  When  it  comes  to  presentation  of 
results,  not  even  the  term  'SIGNAL  TO  NOISE  RATIO*  has  a  universal 
definition.  As  used  here,  signal  to  noise  ratio  is  defined  to  be  ten 
times  the  log  base  10  of  the  ratio  of  the  signal  variance  to  the  error 


variance.  Hie  average  is  again  printed  in  tlie  upper  right-hand  corner 
of  each  plot.  Note  the  broken  T  axis  for  signal  to  noise  ratio. 

Figures  5.9(a)  through  (d)  and  5.10(a)  through  (d)  are  presented  to 
show  the  interaction  of  the  system  gain,  or  signal  to  noise  ratio,  to 
both  the  block  rate  and  data  rate.  For  these  sets  of  plots,  it  would  be 
ideal  if  the  system  signal  to  noise  ratio  vas  not  a  function  of  data 
rate  or  block  rate.  The  ideal  case  would  be  one  for  which  the  output 
signal  to  noise  ratio  would  be  constant  regardless  of  the  data  rate  or 
block  rate.  That  is,  the  quality  of  the  output  images  should  remain 
constant  regardless  of  what  the  input  sequence  contains.  Neither  of  the 
sets  of  plots  show  a  tendency  to  disprove  this  idea. 

The  previous  three  sets  of  plots  were  functions  of  the  overall  system 
gain.  One  of  the  internal  parameters  not  measurable  from  the  output  of 
the  system  that  has  proved  informative,  is  the  predictor  gain.  The 
predictor  gain  is  defined  to  be,  ten  log  base  10  of  the  ratio  of  the 
variance  of  the  signal  to  the  variance  of  the  error  of  prediction.  This 
prediction  error  of  the  system  is  that  error  which  would  be  generated  if 
no  residual  were  available  for  correction.  Recall  that  this  is  one  of 
the  internal  parameters  used  to  judge  what  the  data  rate  of  the  residual 
should  be.  The  plots  of  the  predictor  gain  as  a  function  of  frame 
number  are  provided  in  figure  5.11(a)  through  (d).  Note  that  these 
plots  tend  to  be  much  more  erratic  and  varying  than  do  those  of  the 
entire  system  signal  to  noise  ratio.  Although  not  plotted,  the 
difference  between  system  gain  and  the  predictor  gain  shows  the  gain 
that  has  to  be  generated  by  the  quantizer  and  error  coding  section  of 
the  system. 

Unlike  the  overall  system  gain,  it  is  expected  that  the  predictor 
gain  would  be  related  to  both  the  data  rate  and  block  rate.  Figures 
5.12(a)  through  (d)  show  the  relationship  between  the  predictor  gain  and 
the  block  rate.  As  would  be  expected,  the  points  tend  to  cluster  about 
a  line  with  a  slightly  negative  slope.  Said  another  way,  the  predictor 
gain  is  higher  for  a  smaller  number  of  blocks  in  motion.  The  greater 
the  number  of  blocks  in  motion,  the  lower  the  predictor  gain.  Likewise 
the  predictor  gain  and  system  data  rate  exhibit  similar  characteristics 
in  figures  5.13(a)  through  (d).  The  lower  the  predictor  gain,  the  more 
data  that  has  to  be  transmitted  through  the  quantizer  and  hence,  the 
higher  data  rates. 


Hie  final  plots  provide  s  sort  of  ststistiosl  view  of  the  predictor 
error  signal.  Figures  5.14(a)  through  (d)  present  a  log  histogram  of 
the  prediction  errors.  Figures  5.15(a)  through  (d)  provide  an 
alternative  view  of  the  prediction  error  in  a  log  cumulative 
distribution  function  of  the  error. 


FIGURE  5.5(a)  SL0V  PHONE  1.  Block  Rate 


FIGURE  5.5(b)  SLOW  PHONE  2.  Block  Rate 
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FIGURE  5.5(e)  SLOV  PHONE  3.  Block  Rate 


FIGURE  5.5(d)  FAST  PHONE.  Block  Rate 
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FIGURE  5.7(a)  SLOT  PHONE  1.  Data  Rate  to  Block  Rata 
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FIGURE  5.7(b)  SLOW  PHONE  2.  Data  Rate  to  Block  Rate 
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FIGURE  5.7(c)  SLOW  PHONE  3.  Data  Rate  to  Block  Rate 


FIGURE  5.7(d)  FAST  PHONE.  Date  Rate  to  Block  Rate 
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FIGURE  5.8(a)  SLOW  PHONE  1.  System  Sigaal  To  Noisa  Ratio 
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FIGURE  5.8(b)  SLOF  PHONE  2.  Systea  Signal  To  Noise  Ratio 
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FIGURE  5.8(c)  SLOV  PHONE  3.  System  Signal  To  Noise  Ratio 
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FIGURE  5.8(d)  FAST  PHONE.  System  Signal  To  Noise  Ratio 
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FIGURE  5.9(b)  SLOW  PHONE  2.  Syatea  SNR  to  Block  Rate 


FIGURE  5.9(c)  SLOT  fSGNE  3.  Systea  SNR  to  Block  Rate 


FIGURE  5. 9(d)  FAST  HIGNE.  Systea  SNR  to  Block  Rate 
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FIGURE  5.10(a)  SLOT  PHONE  1.  System  SNR  to  Data  Rate 
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FIGURE  5.10(b)  SLOW  PHONE  2.  System  SNR  to  Data  Sate 
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FIGURE  5.10(c)  SLOF  PHONE  3.  Sy*t*«  SNR  to  Dot*  Rot* 
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FIGURE  5.11(a)  SLOW  PHONE  1.  Predictor  Gaia 
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FIGURE  5.11(b)  SLOW  PHONE  2.  Predictor  Gaia 
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FIGURE  S. 12(a)  SLOV  PHONE  1.  Block  Rata  to  Predictor  Gain 
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The  iteration  methods  ere  also  unable  to  cope  with  ocoluaion  aa  veil  as 
non-translation  types  of  notion.  The  Prediction  Coefficient  method  is 
able  to  partially  alleviate  these  problems. 

The  results  presented  show  a  good  degree  of  data  compression,  while 
at  the  same  time  producing  high  quality  output.  The  algorithm  has  been 
tested  on  data  similar  to  that  which  would  be  encountered  in  a  video¬ 
telephone  setting.  A  large  amount  of  the  algorithm  analysis  has  been 
performed.  To  some  degree  it  has  been  shown  how  various  system 
parameters  effect  the  quality  of  the  output  and  the  data  rate  required 
to  achieve  it.  Also,  the  interaction  among  the  predictor  gain,  the  gain 
produced  by  the  quantizer,  and  the  overall  system  gain  have  been 
investigated.  Actual  output  images  have  been  presented  for  visual 
comparison  of  the  different  preset  parameters.  Ideas  and  suggestions 
were  provided  for  possible  future  research  to  continue  with  the  work 
that  has  been  started  here. 
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requirement  in  background  and  shadow  areas  than  wonld  be  expected.  In 
these  background  and  shadow  areas,  both  the  signal  variance  and  signal 
power  are  low.  There  say  be  a  way  to  take  into  account  both  the  signal 
variance  and  the  signal  power  to  improve  the  quantizing  quality. 

Finally,  work  is  needed  in  the  area  of  noise  effects.  A  study  of  how 
noise  in  the  input  sequence  effects  both  the  prediction  process  and  the 
output  quality  of  the  image  sequence  is  needed.  More  importantly,  how 
noise  in  the  digital  channel  effects  the  operation  of  the  algorithm,  as 
well  as  the  quality  of  the  final  output,  needs  further  research.  The 
implementation  here  assumed  a  noise  free  channel,  but  for  real  world 
applications  an  additive  noise  channel  will  need  to  be  modelled. 

Only  a  few  of  the  many  possible  refinements  and  extensions  of  this 
one  particular  method  for  motion  compensated  image  coding  have  been 
mentioned.  The  entire  field  of  sration  compensated  image  coding  is  still 
in  its  infancy,  with  an  enormous  amount  of  work  remaining.  To  even 
attempt  to  list  all  of  the  possible  paths  for  future  research  would  be  a 
major  undertaking.  The  fields  of  research  remain  wide  open. 

5.3  CONCLUSIONS 

The  method  of  Prediction  Coefficient  Energy  Concentration  has  proved 
to  yield  an  improvement  in  motion  compensated  image  coding  when  compared 
with  previously  published  results.  The  improvements  come  in  both  the 
areas  of  algorithm  complexity  and  in  the  method  for  determining  the 
displacement  estimation.  Previous  methods  for  displacement  estimation 
have  relied  on  an  iterative  procedure  that  was  both  time  consuming  and 
rigidly  fixed  to  translation  types  of  motion  only.  The  method  of 
Prediction  Coefficient  Energy  Concentration  has  replaced  the  iterative 
estimation  procedure  with  a  two-step  estimation  procedure.  The  first 
step  produces  an  estimate  for  the  integer  portion  of  the  displacement. 
The  second  step,  if  needed,  produces  a  set  of  predictor  coefficients 
that  are  able  to  perform  the  same  function  as  the  non- integer  portion  of 
the  displacement. 

The  algorithm  does  not  have  the  problems  of  convergence  that  many  of 
the  previous  methods  exhibit.  The  replscement  of  the  iteration 
procedure  with  the  two-step  prediction  process  relieves  this  problem. 


Further  work  is  required  to  determine  the  effectiveness  of  this 
procedure  ss  veil  ss  how  this  setting  effects  the  actusl  quality  of  the 
output.  Signal  to  noise  ratio  is  a  ten  that  is  easily  calculated,  but 
■ay  not  have  a  direct  relationship  to  the  actual  quality  of  the  output 
iaage  sequence.  Perhaps  some  other  figure  of  merit  would  function 
better. 

The  area  of  the  procedure  that  requires  the  majority  of  the  bandwidth 
is  also  the  area  that  offers  the  greatest  potential  for  further 
research.  This  is  the  area  of  error  quantization.  The  system  as 
presented  here,  assumed  the  error  distribution  to  be  Gaussian  and 
attempted  to  minimize  the  sum  of  the  absolute  value  of  the  error.  This 
is  a  major  oversimplification  of  the  data  that  actually  results.  It  is 
known  that  the  distribution  is  not  Gaussian,  but  estimates  of  the  actual 
distribution  arr  unknown.  Better  estimates  of  the  actual  error 
distribution  are  required.  Perhaps  one  solution  would  be  to  define  a 
set  of  possible  distributions  and  determine  which  of  these  distributions 
the  current  residual  data  best  fits.  It  is  expected  that  the  errors 
that  arise  from  blocks  located  at  the  moving/non-moving  interface  would 
differ  substantially  from  those  blocks  that  tend  to  have  a  constant 
translation.  Yet  another  distribution  would  be  required  for  areas  where 
background  is  uncovered.  If  it  could  be  determined  from  which  area  a 
block  originated,  an  estimate  for  that  local  error  distribution  might  be 
generated. 

Another  problem  in  the  ares  of  quantization  that  merits  further 
research,  is  that  of  possible  data  dependent  thresholding  or 
quantization.  It  has  already  been  stated  that  the  eye  will  allow  a 
higher  degree  of  degradation  for  portions  of  the  image  that  are 
undergoing  translation  than  it  will  for  non-moving  portions.  It  also 
seems  that  the  human  visual  system  is  very  critical  about  small  errors 
when  they  occur  in  areas  of  the  image  with  small  signal  variance.  This 
problem  was  partially  addressed  in  that  the  signal  to  noise  ratio  was 
defined  as  ten  log  the  ratio  of  the  signal  variance  to  the  error 
variance,  as  opposed  to  ten  log  the  ratio  of  the  signal  power  to  the 
error  power  as  is  normally  used.  This  forced  a  better  residual  update 
in  the  areas  of  the  image  with  small  signal  variance,  such  as  the  face 
and  hand  areas  of  the  images  presented.  It  also  forced  a  higher  update 
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alleviated  with  the  help  of  a  low  pass  filter  on  the  output  stage  of  the 
motion  compensation  system  to  smooth  the  output. 

5.2  BBmmfRNTt ATT ONS  FOR  FUTURE  RESEARCH 

This  work  has  only  scratched  the  surface  as  far  as  motion  compensated 
image  coding  is  concerned.  The  groundwork  material  has  been  set  and  a 
path  started,  but  many  sore  questions  have  been  asked  than  have  been 
answered.  Much  of  the  algorithm  analysis  and  implementation  work  yet 
remains. 

One  area  of  the  algorithm  analysis  that  needs  to  be  studied  in  more 
detail,  is  the  area  of  optimal  threshold  and  parameter  selection.  As 
pointed  out  in  the  last  section,  effort  is  needed  to  determine  optimal 
values  of  the  thresholds  and  how  the  choice  of  one  effects  the  output, 
as  well  as  the  action  of  some  of  the  other  parameters.  Ideally,  some  of 
the  parameters  should  be  tied  to  the  quality  of  the  image  data  itself. 
It  will  prove  unproductive  to  try  to  generate  an  output  image  sequence 
with  a  signal  to  noise  ratio  greater  than  that  of  the  input  sequence. 
Along  the  same  lines,  the  difference  in  quality  of  display  devices  can 
also  contribute  to  the  setting  of  various  system  parameters.  That  is, 
the  dynamic  range  and  distortion  of  the  display  device  needs  to  be 
considered  when  the  image  quality  requirements  are  defined.  In  general, 
the  choice  for  the  threshold  for  determining  if  the  displacement  is 
integer  or  non-integer  as  well  as  the  choice  for  the  threshold  for 
determining  if  sub-block  processing  is  required,  requires  some  more 
effort.  These  two  thresholds  are  not  unrelated.  For  example,  if  the 
error  threshold  is  increased  for  the  integer/ non-integer  test,  allowing 
more  blocks  to  pass  as  integer  displacement,  fewer  will  even  reach  the 
portion  of  the  algorithm  that  tests  for  sub-block  processing.  Some  work 
needs  to  be  performed  to  determine  what  the  optimal  threshold  settings 
are,  how  they  are  related  to  the  actual  image  data,  and  also  how  they 
are  related  to  one  another. 

One  of  the  other  preset  system  parameters  that  warrants  further 
research,  is  the  value  for  the  output  signal  to  noise  ratio.  As  was 
noted  in  one  of  the  previous  chapters,  when  a  fixed  length  buffer  is 
used  in  the  implementation,  the  value  for  this  required  signal  to  noise 
ratio  of  the  output  could  be  used  as  a  buffer  regulating  parameter. 


5.1  aBMATWTWn  PBOWTEMS 


The  overall  results  of  this  method  for  motion  compensated  image 
coding  have  proved  promising.  A  fairly  high  quality  output  has  been 
achieved  at  a  nominal  overall  data  rate.  Still,  for  video-telephone 
applications  the  data  rate  remains  too  high.  There  remains  a  number  of 
problems  with  both  the  algorithm  itself  and  its  implementation. 

Perhaps  the  biggest  problem,  other  than  the  required  bandvidth,  for  a 
video-telephone  implementation,  is  that  of  computational  speed.  It  is, 
and  would  be,  impractical  to  try  to  implement  the  algorithm  in  anything 
other  than  a  parallel  architecture  hardware  device.  At  various  points 
throughout  the  different  chapters,  sections  of  the  algorithm  that  lend 
themselves  to  parallel  implementation  were  noted.  Real  time  processing 
of  video  images  is  a  very  calculation  intensive  undertaking  that 
interfaces  well  with  parallel  processing  techniques.  The  problem  of 
processing  time  may  be  negligible  when  a  real  time  hardware  parallel 
processor  is  used. 

Another  problem  of  the  computer  simulation  was  that  of  the  choice  of 
the  various  thresholds  and  other  system  parameters.  No  attempt  was  made 
to  optimize  all  of  these  parameters.  Some  of  the  thresholds  and 
parameters  were  chosen  to  estimate  the  noise  introduced  into  the  system 
by  the  imaging  apparatus,  while  others,  such  as  the  predictor 
coefficient  coding  length,  were  chosen  by  a  mathematical  model.  The 
image  quality  and  data  rate  requirements  are  opposing  ends.  That  is, 
minimizing  the  data  rate  tends  to  decrease  the  image  quality,  as  does 
increasing  the  image  quality  require  an  increase  in  data  rate.  Both  of 
these  are  very  dependent  upon  the  settings  of  the  various  thresholds. 

One  of  the  problems  that  remains  in  the  output  image  can  only  be  seen 
when  the  images  are  viewed  in  sequence  as  they  were  intended  to  be. 
This  problem  comes  about  because  of  the  ability  of  the  human  eye  and 
visual  system  to  detect  very  minute  changes  in  gray  level,  if  it  has 
some  geometric  structure.  This  can  be  seen  to  occur  in  the  image 
sequence  at  lower  data  rates  when  the  outlines  of  the  sub-blocks  of  the 
images  become  visible.  Because  of  the  straight  line  geometric  nature  of 
the  block  boundaries,  they  appear  more  noticeable  than  one  would  expect 
for  the  size  of  the  actual  error.  Perhaps  this  could  be  partially 
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Appendix  A 

PEL  RECURSIVE,  COEFFICIENT  RECURSIVE  AND  RESIDUAL 
RECURSIVE  DISPLACEMENT  ESTIMATION 


DISPEST1 


PROGRAM  DISPEST1 


PROGRAM  DISPEST1  AUTHOR  -  CARL  BOWLING  DATE  9/12/83 

PROGRAM  FOR  DISPLACEMENT  ESTIMATION 

METHOD  IS  COEFFICIENT  RECURSIVE  BY  STULLER  AND 
NETRAVALI  OF  BELL  LABORATORIES.  IT  MAY  ALSO  BE  USED 
FOR  PEL  RECURSIVE  BY  NETRAVALI  AND  ROBBINS 

SUBROUTINES  OR  FUNCTIONS  RBQUIRED  - 

XFORM  -  ROUTINE  TO  PERFORM  TRANSFORMATION 
PLOTS  -  CAL COMP  TRANSLATOR  AND  PLOT  SOFTWARE 
RINTRP  -  FUNCTION  FOR  INTEGER  INTERPOLATION 


REAL  V(2,2)  ,H(  8,8)  ,CLL<6,8)  ,CP(6.256)  ,CPP(6,8)  ,CPPP(6,8) 
INTEGER  IL(6,256) ,IP(6,256) ,ILL(6,8) ,IPP(6,8) 

REAL  G( 6,8) ,GG(6,8) ,CLLL(6.8) ,ERRQR(6,8) ,DISP(250) ,XIT(250) 
DATA  PI/3.14159265/ 

INITIALIZE  MATRICES  TO  VALUE  OF  128 

DATA  JL! 1536*128/, IP/1536*128/ 

SET  UP  THE  TRANSFORMATION  MATRIX  IN  THE  H  ARRAY 


DATA  H/ 12*1. 00, 4*— 1 . 00, 2*1. 00, 4*— 1 .00,4*1.00,2*-1 .00, 
2*1. 00, 2*— 1.00,  1.00.2*-1.00,2*1.00,2*-1.00, 
2*1.00,2*-1 .00,  1.00,  -1.00,2*1.00,  -1.00, 
1.00,  -1.00,  1.00,2*-1.00,  1.00,  -1.00, 
2*1.00,  -1.00,  1.00,  -1.00,  1.00,  -1.00, 
1.00,  -1.00/ 


METHOD  IS  USED  TO  DETERMINE  IP  PEL  RECURSIVE  OR  COEFFICIENT 
RECURSIVE  DISPLACEMENT  ESTIMATION  IS  USED.  APPROPRIATE  LINES 
MUST  ALSO  BE  COMMENTED  OUT  DUE  TO  THE  COMPILER  USED. 

METHOD  -  0  — >  PEL  RECURSIVE 
METHOD  -  1  — >  COEFFICIENT  RECURSIVE 

LOOP  IS  USED  TO  DETERMINE  IF  EVERY  LOOP  ITERATION  OR 
EVERY  BLOCK  ITERATION  IS  TO  BE  USED. 

LOOP  -  0  — >  EVERY  BLOCK  ITERATION 
LOOP  -  1  — >  EVERY  LOOP  ITERATION 

NPLOT  IS  A  FLAG  TO  DIRECT  THE  PLOTTING  OF  THE  DATA 
NPLOT  -  0  — >  GO  TO  4025  SCREEN 
NPLOT  «  1  — >  GO  ID  4662  PLOTTER 

LOOP  -  1 
METHOD  -  1 
NPLOT  *  0 

SQRT2  -  l./((2.)**.5) 

IF ( METHOD. BQ.O) SQRT2  -  0.0 
V(l.l)  =«  SQRT2 
V(1.U  -  V(l,l) 

V(2.1)  =  V(l.l) 

V(2,2)  -  -V(l,l) 

FAC  -  . 5  *SQRT2 
DO  10  1-1,8 
DO  10  J-1,8 
H(I, J)  -  H(I, J)*FAC 
10  CONTINUE 

IF ( METHOD. NE.O)  GO  TO  40 
DO  20  >1,8 
20  H( I, I)  -  1.0 
DO  30  >1,2 
30  V(I,I)  -  1.0 
40  CONTINUE 


C************* ****••*•* ********** •**••*****•* ****** *••*•*••* 

C*  GENERATE  THE  TEST  ARRAYS  IN  IP  AND  EL  * 

C*  IP  IS  THE  PRESENT  IMAGE  FRAME  * 

C*  IL  IS  THE  PREVIOUS  IMAGE  FRAME  * 

C*  THE  TEST  ARRAY  IS  A  RADIALLY  DECAYING  FREQUENCY  • 

C*  MODULATED  COSINE  FUNCTION.  * 

Cl*********************************************************** 

C 

DO  50  1-1.6 
DO  50  J*1 ,60 

R-SQRT(  FLOAT (  (1+7)  **2  +  J*J  )  ) 

IF  (  R  .GT.  60.0  )  GO  TO  50 
P-(l.  -  R/60.)*10.  +  10. 

IPT-100.  •  EXP(-.01*R)  •  COS(2.*PI*R/P)  +  128. 
IP(I,J+128)«IPT 
IP(I,129-J) -IPT 
IL(I, J+126)-IPT 
IL(I, 127-J) «IPT 
50  CONTINUE 
C 

C*********************************************************** 
C*  * 

C*  TRANSFORM  THE  ORIGINAL  ARRAY  WITH  SUBROUTINE  XFORM  * 
C*  • 

C********* *•*•****•• *********************** ***************** 

C  NR  -  THE  NUWER  OF  ROWS  IN  TRANSFORMATION 

C  NC  -  THE  NUWER  OR  COLUMNS  IN  TRANSFORMATION 

NR-6 
NC-8 
C 

C  TRANSFORM  IS  (V)*(I)*(H) 

C 

DO  80  J-1,256.8 
C 

C  J  IS  THE  BLOCK  NUWER 

C 

DO  60  II  -  1.6 

C  U  IS  THE  ROW  NUWER 

J7  -  J  +  7 
13-0 

DO  60  12  -  J.J7 
13  -  13  +  1 

CPPP(K1.I3)  -  IP(I1,I2) 

60  CONTINUE 

CALL  XFORM(CPP.V.CPPP,H.NR,NC) 

DO  70  II  -  1,6 
13  -  0 

DO  70  K2-J.J7 

13  -  13  +  1 

CP(K1,I2)  -  CPPdl  ,13) 

70  CONTINUE 
80  CONTINUE 
C 

C*********************************************************** 

c*  • 

C*  END  TRANSFORMATION  -  START  OF  DISPLACEMENT  ITERATION* 

C*  * 

C*********************************************************** 

c 

C  EPSIL  -  GAIN  FACTOR  (READ  IN  INTERACTIVELY) 

C  DA  IS  THE  ACTUAL  FRAME  DISPLACEMENT  IN  PIXELS 

C  DI  IS  THE  ESTIMATE  OF  THE  PIXEL  DISPLACEMENT 

C  ITNUM  IS  THE  ITERATION  LOOP  COUNTER 

C 

DA  -  2.0 
XOFF  -  1.35 
YOFF  -  1.5 

IF(NPLOT.  BQ.O)  CALL  PLOTSdBUF.l  ,15) 

IF ( NFLOT. BQ .0 ) XOFF  -  .1 
IF(NFLOT. BQ.O) YOFF  -  .1 
90  CONTINUE 
WRITE( 6,500) 

READ (5, 510) EPSIL 
IF(EPSIL.GT.l)  STOP 
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IF(NPLOT. EQ.O) CALL  ERASE 
DI  *  0.0 

I  IS  THE  BLOCK  NUMBER  (USE  ONLY  10TH  -  23RD) 

ITNUM  =  0 
DO  150  1*10,23 

JJ  -  START  COLUMN  LOCATION  OF  CURRENT  BLOCK 
JD  -  START  COLUMN  LOCATION  OF  PREVIOUS  BLOCK 

JJ  -  ( 1-1) *8  +  1 
J  *  1 
K  *  3 

IF ( LOOP . BO . 0 )  GO  TO  95 

IF  LOOP  SET  TO  0,  COMMENT  OUT  THE  FOLLOWING  TWO  LINES 
SOME  COMPILERS  MAY  NOT  ALLOW  THE  JUMP  AROUND  DO  LOOPS 

DO  140  J-1,8 
DO  140  K*3 ,4 
RD  -  FLOAT (JJ)  -  DI 
JD  -  RD 

RDIFF  -  RD  -  FLOAT (JD) 


CALCULATE  INTERPOLATED  DISPLACED  BLOCK 


JD7  -  JD*7 
JDD  *  0 

DO  100  L-JD.JD7 
JDD  -  JDD  +  1 
DO  100  M-1.6 

X  -  RI(IL<M.L),IL<M.L+1). RDIFF) 
CLLL(M.JDD)  *  X 


END  OF  FRAME  INTERPOLATION  DETERMINATION 
FIND  THE  FRAME  ERROR 


CALL  XFORM( CLL, V,  CLLL,  H, NR, NC) 

DO  110  L*l,8 
DO  110  1H3  „4 

ERROR(M.L)  -  CP(M, L+JJ-1)  -  CLL(M.L) 
CONTINUE 


FIND  TIME  DOMAIN  GRADIENT  OF  THE  DISPLACED  FRAME 
CLLL.  GRADIENT  FOUND  BY  CENTRAL  DIFFERENCE  METHOD 


DO  130  L-3.4 
DO  120  M-2,7 

GRAD1  *  CLLL(L+1,M)  -  CLLL(L-l.M) 
GRAD1  -  0. 

GRAD2  -  CLLL(L, M+l)  -  CLLL(L.M-l) 
G(L,M)  *  ( GRAD1  +  GRAD2)  *  .5 
GRAD1  -  <CLLL(L,2)  -  CLLL(L, 1) ) *2. 
GRAD2  -  <CLLL(L+1,1)  -  CLLL (L- 1,1) ) 

GRAD2  -  0. 

G(L, 1)  -  (GRAD1  +  GRAD2)  *.5 
GRAD1  -  (CLLL(L,8)  -  CLLL(L, 7) ) *2 . 
GRAD2  -  (CLLL(L+1,8)  -  CLLL(L-1,8) ) 
GRAD2  -  0. 

G(L,8)  *  (GRAD1  +  GRAD2)*.5 
CONTINUE 

CALL  XFORM(GG,V,G,H,6 ,8) 
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c* 

C*  GENERATE  THE  DISPLACEMENT  ESTIMATE 

C* 


IF(LOOP.BQ.l)  GO  TO  135 
C 

C  IF  LOOP  IS  SET  TO  1,  COMMENT  ODT  THE  FOLLOWING  TWO  LINES 

C  SOME  COMPILERS  MAT  NOT  ALLOW  THE  JUMP  AROUND  DO  LOOPS 

C 

C  DO  140  J  -  1,8 

C  DO  140  K  =»  3.4 

135  DI  -  DI  -  EPSIL*ERROR(K, J) *GG(K, J) 

ITNUM  »  ITNUM  +  1 
XIT(ITNUM)  *  ITNUM 
DISP( ITNUM)  -  DA  -  DI 
140  CONTINUE 
150  CONTINUE 

IF ( NFLOT. EQ.l) CALL  PLOTS  ( IB UF, 1,1 5) 

CALL  PLOTtXOFF, Y0FF.-3 ) 

IF ( NPLOT. EQ . 0 ) XOFF  =0.0 
IF ( NPLOT. BQ . 0 ) TOFF  =  0.0 
CALL  FACTOR* .65) 

IF ( NPLOT.  EQ.O) CALL  FACTOR* .25) 

CALL  SCALE*DISP, 5 . ,80 ,1) 

CALL  SCALE(XIT,8 . ,80,1) 

CALL  AXIS (0.0, 0.0, 'ITERATION  NUMBER' ,-16 .8 . ,0 .0 ,XIT(81) ,XIT( 82) ) 
CALL  AXIS (0..0., 'DISPLACEMENT  ERROR' ,18,5 . ,90 . ,DISP( 81) ,DISP( 82) ) 
IF ( NPLOT. B8.0)  GO  TO  155 
C  IF (NPLOT. EQ.l) CALL  NEWPEN* 'BLAC' ) 

WRITE* 6, 520) 

520  FORMAT* IX, '/•  HIT  ANT  SINGLE  DIGIT  NUMBER  AND  RETURN  IF  READY' ) 
READ* 5. 530) IANS 
530  P0RMATU2) 

C  IF ( NPLOT. BQ . 1 ) CALL  NEWPEN* 'BLAC' ) 

155  CONTINUE 

CALL  LINE(XIT, DISP, 80,1,1,3) 

CALL  SYMBOL* 6. 3, 4. 45,. 20, 44, 0..-1) 

CALL  STMBOL(6. 5. 4. 5,. 14. '-',0.,1) 

CALL  NUMBERS. 7,4  .5.. 14. EPSIL, 0.0, 5) 

IF ( LOOP . BQ . 1 ) CALL  SYMBOL(2. ,5 . , .2 . ' EVERY  LOOP  ITERATION  '.0..21) 
IF ( LOOP . BQ . 0 ) CALL  SYMBOL* 2. ,5 ., .2, 'EVERY  BLOCK  ITERATION' ,0 . ,21) 
IF(METHOD. EQ.l) CALL  SYMBOL(2. ,5 .3 , .2 , ' COEFFICIENT  RECURSIVE' ,0 . 

#  .21) 

IF ( METHOD. BQ.O) CALL  SYMBOL(2. ,5 .3 . .2 , ' PEL  RECURSIVE' ,0 . ,13) 

IF(NFLOT.BQ.l)  GO  TO  90 

CALL  TSEND 

CALL  AN MODE 

GO  TO  90 

160  CALL  PLOT*  10.. 10., 999) 

500  FORMAT* IX,  7*'.  'WHAT  VALUE  FOR  THE  GAIN  EPSILON?  ') 

510  F0RMAT(F10.6) 

STOP 

END 
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PROGRAM  RESDISP 

PRT  FILE  5279  FROM  VS1  COPT  001  NOHOLD 

C*********************************************************** 

c*  * 

C*  PROGRAM  BY  -  CARL  BOWLING  * 

C*  LATEST  UPDATE  -  OCTOBER  26,  1983  * 

C*  * 

C*  PROGRAM  FOR  DISPLACEMENT  ESTIMATION  * 

C*  METHOD  IS  RESIDUAL  RECURSIVE  * 

C*  * 

C*  SUBROUTINES  OR  FUNCTIONS  REQUIRED  -  * 

C*  * 

C*  XFORMB  -  PERFORMS  A  1  OR  2  DIMENSIONAL  UNITARY  * 

C*  TRANSFORMATION  ON  THE  INPUT  ARRAY.  * 

C*  COMP  -  ADAPTIVE  HYBRID  PICTURE  CODING  DATA  * 

C*  COMPRESSION  ROUTINE.  * 

C*  PLOTS  -  ALL  THE  CAL  COMP  DRIVER  SOFTWARE  AND  THE  * 

C*  TRANSLATOR  SOFTWARE  FOR  THE  TEKTRONIX  * 

4662  PLOTTER  AND  4025  DISPLAY.  • 

C*  RINTRP  -  FUNCTION  TO  INTERPOLATE  BETWEEN  ADJACENT  * 

C*  INTEGER  VALUES.  * 

C*  * 

(;•*••••*••*«••*«*•*••*****••«*••••**••**•**•*•*••*••******•* 

C  _ _ 

EXTERNAL  HADGEN 

REAL  DLIF( 8,8) ,DLTF(8,8) , ERROR( 8.8) , CP( 8.256) 

REAL  COEFF( 8.3) ,CL( 8.256) ,XIT(250) ,DISP(250) 

INTEGER  IL(8,256) ,IP(8,256) 

COMMON  / CODPRM/  IR. IC,  EBSZ.NTYPE.CBSZ.OTYPE. IBLKSZ 
COMMON  /FT/  KFFT 

DATA  PI/ 3. 141 5926 5/  _ 

DATA  NTYPE1/ 'IN*4 '/ .OTYPE1/ 'REAL'/ 

C 

C  NTYPE1  IS  THE  DATA  TYPE  OF  THE  INPUT  ARRAY  (INTEGER* 4) 

C  OTYPE1  IS  THE  DATA  TYPE  OF  THE  OUTPUT  ARRAY  (REAL*4) 

C 

C  INITIALIZE  PAST  AND  PRESENT  FRAME  MATRICES  TO  VALUE  OF  128 

C 

DATA  IL/ 2048* 128/. IP/ 2048* 128/ 

C 

C  SET  CONSTANTS  FOR  TRANSFORM 

C 

NTYPE  «  NTYPE1 
OTYPE  -  0TYPE1 
KFFT  -  0 
IBLKSZ-  8 
13  SZ  -  8 
OBSZ  -  8 
C 

C«M*M»MM«*M***MM«***MMM*«M*M*MMM»*M*MM*M* 

c*  • 

C*  GENERATE  THE  TEST  ARRAYS  IN  IP  AND  IL  * 

C*  IP  IS  THE  PRESENT  IMAGE  FRAME  * 

C*  IL  IS  THE  PREVIOUS  IMAGE  FRAME  * 

C*  THE  TEST  ARRAY  IS  A  RADIALLY  DECAYING  FREQUENCY  * 

C*  MODULATED  COSINE  FUNCTION.  ONLY  A  PORTION  OF  THE  * 

C*  ENTIRE  IMAGE  IS  USED.  * 

C*  * 

C******************************************** *************** 

c 

DO  10  1-1,8 
DO  10  J-1,60 

R-SQRT(  FLOAT (  (1+7) **2  +  J*J  )  ) 

IF  (  R  .GT.  60.0  )  GO  TO  10 
P-(l.  -  R/ 60. ) *10.  +  10. 

IPT-100 . *EXP(-.01*R) *COS( 2 . *PI*R/P)  +  128. 

IP(I, J+128)-IPT 
IP( 1,129-J) -IPT 
IL(I, J+126)-IPT 
IL(I,127-J) -IPT 
10  CONTINUE 
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C*****«****»****»«*M**M**«**M**««*»*MM**M*M*MM***** 

c*  * 

C*  START  OF  DISPLACEMENT  ITERATION  * 

C*  • 

C*********************************** *************** ********* 

C 

C  EPSIL  -  GAIN  FACTOR  (READ  IN  INTERACTIVELY) 

C  DA  IS  THE  ACTUAL  FRAME  DISPLACEMENT  IN  PIXELS 

C  DI  IS  THE  ESTIMATE  OF  THE  PIXEL  DISPLACEMENT 

C  ITNUM  IS  THE  ITERATION  LOOP  COUNTER 

C 

DA  -  2.0 
20  WRITE( 6,30) 

30  FORMAT( IX, '/*', 'WHAT  VALUE  FOR  EPSILON?  F10.6  FORMAT') 
READ (5, 40) EPSIL 
40  FORMAT(F10.6) 

C 

C  CHECK  IF  GAIN  IS  GREATER  THAN  1.  IF  SO  STOP  PROGRAM 

C 

IF(EPSIL. GT.l .)  GO  TO  130 
CALL  PLOTS (IBUF, 1,9) 

DI  =  0.0 
C 

C*************************************** ******************** 

c*  * 

C*  I  THE  BLOCK  NUMBER  (USE  ONLY  10TH  -  23RD)  * 

C*  ITNUM  -  THE  ITERATION  NUMBER  * 

C*  JJ  -  START  COLUMN  LOCATION  OF  CURRENT  BLOCK  * 

C*  JD  -  START  COLUMN  LOCATION  OF  PREVIOUS  BLOCK  * 

C*  * 

C*********** **•***•***••*•*•*•* *********************** ****** 

c 

ITNUM  -  1 
XIT(l)  -  1. 

DISP(l)  -  2.0 
DO  120  1-10,23 
JJ  -  ( 1-1) *8  +  1 
DO  120  J  -  1.8 
DO  120  K  -  3.3 
RD  -  FLOAT  (JJ)  -  DI 
JD  -  R D 

RDIFF  -  RD  -  FLOAT(JD) 

C 

C*********************************************************** 

c*  * 

C*  CALCULATE  INTERPOLATED  DISPLACED  BLOCK  * 

C*  IN  THE  TIME  DOMAIN  * 

C*  * 

C*************************************************** ******** 

c 

JD7  -  JD+7 
JDD  -  0 

DO  50  L-JD.JD7 
JDD  «  JDD  +  1 
DO  50  M-1,8 

X  -  RINTRP( IL(M, L) ,  IL(M, L+l) , RDIFF) 

50  DLIF(M.JDD)  =*  X 
DO  60  L-1,8 

60  WRITE(l) (DLIF(L.M) . M-1,8) 

REWIND  1 
C 

C*********************************************************** 

c*  * 

C*  END  OF  DISPLACED  BLOCK  INTERPOLATION  * 

C*  * 

C*  CALCULATE  THE  RESIDUAL  SEQUENCE  FOR  USE  AS  THE  * 

C*  GRADIENT.  * 

C*  * 

C************************************************** ********* 

c 

NTYPE  »  0TYPE1 
IR  -  0 
IC  ■  +1 

CALL  XFORMB(HADGEN, 8 ,8 , 1 ,2 , DLTF) 

REWIND  2 
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DO  70  L=l,8 

70  READ(2) (DLTF(L.M) , M=1 ,8) 

REWIND  2 

CALL  COMP(DLTF,8 ,8,3 ,COEFF) 

REWIND  1 
DO  80  L=1 ,8 

80  WRTTE(l) (DLTF(L.M) ,Mhl,8) 

REWIND  1 
IR  =  0 
IC  =  -1 

CALL  IFOEMB ( HADGEN , 8 , 8 , 1 , 2 . DLTF ) 

REWIND  1 
REWIND  2 
DO  90  L*l,8 

90  READ(2) (DLTF(L.M) ,M=1,8) 

DO  100  L=1 ,8 
DO  100  M-1,8 

100  DLTF(L.M)  =  DLTF(L.M) / FLOAT ( IBLKSZ) 

REWIND  1 
REWIND  2 
C 

(:«*•*••••*****•*••*•*******•«•**•**•••*•*••**•**•****•**•••* 

c*  * 

C*  END  OF  RESIDUAL  CALCULATION  * 

C*  * 

C*  FIND  THE  TIME  DOMAIN  FRAME  ERROR  • 

C*  * 

C*********************************************************** 

C 

DO  110  L=*l  ,8 
DO  110  M-3.3 

110  ERROR(M.L)  =  FLOAT( IP(M, L+JJ-1) )  -  DLIF(M.L) 

C 

C*»oo«mmm«*«**m*mm*m*mm«***m**mmmmm**mmm 

c •  * 

C*  GENERATE  THE  DISPLACEMENT  ESTIMATE  * 

C*  * 

c*mm*m**m»m*«****mmmm«*mmmm*«*m***m**»*m»*** 

c 

DI  -  DI  -  EPSIL*ERROR(K, J) *DLTF(K» J) 

ITNUM  -  ITNUM  +  1 
XIT(ITNUM)  -  ITNUM 
DISP( ITNUM)  *  2.0  -  DI 
120  CONTINUE 
C 

C  ALL  OF  THE  CALLS  BELOW  ARE  USED  TO  PLOT  THE  RESULTS 

C 

CALL  PLOT(1.35,6 .0,-3) 

CALL  FACTOR ( .65) 

CALL  SCALE(DISP. 5 . ,80,1) 

CALL  SCALE( XIT, 8 . , 80 , 1 ) 

CALL  AXIS (0.0. 0.0. 'ITERATION  NUMBER' ,-16 .8 . .0 .0 ,XIT( 81) .XIT( 82) ) 
CALL  AXIS(0. .0. . 'DISPLACEMENT  ERROR' ,18,5 . ,90 . ,DISP( 81 ) ,DISP( 82) ) 
CALL  LINE( XIT, DISP, 80,1,1,3) 

CALL  SYMBOL(6.3,4.45,.20,44,0.,-1) 

CALL  SYMBOL ( 6. 5, 4. 5, .14, '*',0. .1) 

CALL  NUMBER(6 . 7,4.5, .14, EPSIL, 0.0,5) 

CALL  SYMBOLU.  0,5.0,  .28.  'DISPLACEMENT  ESTIMATION'  ,0 .  ,23) 

CALL  SYMBOLU. 0,5. 4,  .20,  'RESIDUAL  RECURSIVE' ,0  .  ,18) 

CALL  PLOT( 10. ,10. ,2) 

CALL  TSEND 
CALL  AN MODE 
GO  TO  20 

130  CALL  PLOT(10. ,10. ,999) 

STOP 

END 
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c 

IERSDM  =  0 
IERSQR  -  0 
K  -  ISTY  -  1 
DO  250  1-1.16 
K  -  K  +  1 
L  =  II  -  1 
DO  250  J=1 ,16 
L  -  L  +  1 

IER  =*  B(K,L)  -  ERROR(I.J) 

EST(K.L)  -  ERROR(I.J) 

ERRSIG(K.L)  -  EB16(I, J) 

IERSDM  -  IERSDM  +  IER 
IERSQR  -  IERSQR  +  IER'IER 
IER  -  ER2,6(I»  J)  +  128 
IF(IER. GT.255) IER-255 
IF(IER.LT.l)  IER  -  1 
HIST(IER)  -  HIST (IER)  +  1. 

250  CONTINUE 
C 

C  CALCULATE  THE  MEAN  AND  VARIANCE  OF  THE  QUANTIZED  VERSION  FOR 

C  COMPARISON. 

C 

A1  -  IERSDM 

SI  -  SQRT( (FLOAT (IERSQR)  -  (A1*A1/ 256 . ) ) /255 . ) 

A1  -  Al/256. 

IF(IPRT. BQ. 2) VRITE( 6,1060) IBLK, JBLK, MIN, ( LOC(KC) ,KC=1,2) .AVG.STD, 
#  Al.Sl.INTBIT.WrrS 
260  CONTINUE 
C 

C  CALCULATE  THE  SIGNAL  TO  NOISE  RATIO 

C 

ITOTAL  -  0 
DO  270  I-l.ISIZE 
DO  270  J-l. JSIZE 
ITOTAL  -  ITOTAL  +  B(I.J) 

270  CONTINUE 

TOTAL  -  FLOAT (ITOTAL) 

XNEAN  *  TOTAL/ SIZE 

SSE  -  0.0 

SSS  -  0.0 

SSP  -  0.0 

DO  290  J-l, JSIZE 

PSSS  -0.0 

PSSE  -  0.0 

PSSP  -  0.0 

DO  280  I-l.ISIZE 

PSSP  -  PSSP  +  ERRSIGd,  J) *ERRSIG(I, J) 

IERR  -  B(I, J)  -  ESTd.J) 

INTER  -  B(I,J) 

SIG  -  FLOAT  (INTER)  -  MEAN 

PSSS  -  PSSS  +  SIG*SIG 

PSSE  -  PSSE  +  FLOAT(IERR*IERR) 

280  IAd.J)  -  ESTd.J) 

SSS  «  SSS  +  PSSS 

SSP  -  SSP  +  PSSP 

290  SSE  -  SSE  +  PSSE 

IF(IPRT.BQ.2)WRITE(6,1070)SSE, SSS 
SNR  -  10. *ALOG10( SSS/ SSE) 

SNRP  -  10. *AL0G10( SSS/ SSP) 

SNRARK  IMAGE)  -  SNR 
SNRPAR( IMAGE)  »  SNRP 
IF(IPRT.BQ.l) WRITE (6, 1080) SNR, SNRP 
C  _ _ 

C  WRITE  OUTPUT  IMAGE  OUT  TO  UNIT  10 

C 

CALL  COMTAL(EST, ISIZE, JSIZE, 10) 

C  WRITE(6,1090)FF 

IF (I PRT . GE . 1 ) WRITE ( 6 , 1 1 00 ) 

C 

C  CALCULATE  MOTION  RATE  (#B LOCKS  IN  MOTION) 

C 


NTAB  -  0 
DO  330  1-1,24 
DO  325  K-1,16 
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C*  NON- INTEGER  DISPLACEMENT  ASSUMED  AT  THIS  POINT  * 

C*  * 

C*  LOCATE  THE  MINIMUM  QUADRANT  WITH  RESPECT  TO  THE  * 

C*  INTEGER  DISPLACEMENT.  * 

C*  * 

CM******************************************************** 

c 

CALL  FILL( IBLK, JBLK, I1M1, J1M1 ,8 ,MTAB,24,16 , 'RR' ) 

INTBIT  =  1 

CALL  LOCATE(LOC, F,  21 .21 , MIN) 

NBITS  -  WITS  +  21  +  4*NPB ITS 
C 

C  SET  UP  THE  REGRESSION  PROBLEM 

C 

X  »  0 

DO  190  L-1.8 
I  -  L  +  IXOFF 
DO  190  M-1.8 
J  -  M  +  IYOFF 
I  *  X  +  1 
DO  180  N-1.4 

X(K,N)  -  D8( J+QOFF(MIN, N,l) , I  +  QOFF(MIN, N, 2) ) 

180  CONTINUE 

OUT(X)  *  C8(M, L) 

190  CONTINUE 
C 

C - 

c 

C  CALCULATE  THE  NEW  PREDICTION  COEFFICIENTS 

C 

C - 

c 

CALL  COEFGNd,  OUT,  RCOEF,  X) 

C 

C  QUANTIZE  THE  PREDICTOR  COEFFICIENTS 

C 

DO  200  1-1,4 

200  RCOEF(I)  -  (AINT(RCOEFd)  ‘RLEVEL  +  .5))/RLEVEL 
IF(X.NE.129)  GO  TO  210 
RCOEF(l)  -  0.0 
RC0EF(2)  -  0.0 
RCOEF  ( 3 )  -1.0 
RCOEF (A)  -  0.0 
MIN  -  1 
210  CONTINUE 
C 

C*M*MM**MM**MMMM****M*M««****MMMMM**MMM 

c*  * 

C*  CALCULATE  THE  PREDICTED  IMAGE  BLOCK  * 

C*  * 

£•**•**••*«••••***••***•*••*•••**•••**•*•*••**•*•*•******* 

C 

CALL  PRED ( C8 . D8 , E8 , 2  8 , 8 , IXOFF , I YOFF , VA . AV , QOFF , MIN . RCOEF , INTB  IT , 

#  NBITS, RETCOD, ER8 . DIST, LVLARA) 

IFdPRT.  EQ.2)  WRITE(6,1050) II , J1 ,MIN. LOC(l) ,L0C(2) , AV, VA. INTBIT, 

#  NBITS, RETCOD 

IF ( RETCOD. BQ.l) NBITS  =■  WITS  -  21 
IF ( RETCOD. EQ.l)  GO  TO  170 
C 

C  WRITE  ESTIMATE  AND  ERROR  OUT  TO  16X16  ARRAY 

C 

X  =■  IBLK1 
DO  220  1*1,8 
X  *  X  +  1 
L  *  JBLK1 
DO  220  J-1,8 
L  *  L  +  1 

ERROR(K.L)  -  E8d,J) 

ER16 (X,L)  -  ER8(I,J) 

220  CONTINUE 
230  CONTINUE 
240  CONTINUE 
C 

C  WRITE  QUANTIZED  PREDICTION  VALUE  OUT  TO  ESTIMATED  MATRIX 

C  CALCULATE  THE  MEAN  AND  VARIANCE  OF  THE  PREDICTION  ERROR 


146 


C8(K,L)  =■  C(I,J) 

140  CONTINUE 

IAVG  -  IFIX ( FLOAT (ITOTAL)/ 64.  +  .5) 

CALL  INITKD8  .28,28,  IAVG) 

150  CONTINUE 
C 

C  TRANSFER  DATA  FROM  D  ->  D8 

C 

ISTX  *  IX  -  10  +  JBLK1 
ISPX  *  ISTX  +  27 
ITT  *  1ST!  -  10  +  EBLK1 
IPY  ■  ITY  +  27 
ISX  -  ISTX 

mx  -  ispx 

1ST  -  ITY 
IHY  -  IPY 

IF (IB LOCK. BO. 0)  GO  TO  155 
IF(ISX.LT.l)  ISX  -  1 
IF(IHX.GT. ISIZE)  IHX  -  ISIZE 
IF(ISY.LT.l)  ISY  -  1 
IF ( IHY. GT. ISIZE)  IHY  -  ISIZE 
155  M  *  0 
N  -  0 

IF(ISTX. LT.l)  M  *  1  -  ISTX 
IF(ITY.LT.l)  N  -  1  -  ITY 
KCON8(l)  «  N  +  1 
KCON8 ( 2)  -  M  +  1 

K  ■  M 

DO  160  J-ISX.IHX 
K  -  K  +  1 
L  -  N 

DO  160  I- ISY, IHY 
L  -  L  +  1 
D8(L,K)  -  IAd.J) 

160  CONTINUE 
C 

C - 

c 

C  CALCULATE  THE  METRIC  FOR  THE  8  BY  8  BLOCKS 

C 

C - 

c 

C  IF  (IB  LOCK.  HQ  1)  CALL  INITKF,  21 ,21 ,0) 

CALL  METRIC (CS . D8 ,F, 28 .21 ,8 , IMIN, LOC,  IB  LOCK, KCON8 , ITHRS8 
#  11.11) 

IIOFF  -  LOC( 2)  -  1 
IYOFF  -  LOC(l)  -  1 
C 

C  TEST  FOR  INTEGER  DISPLACEMENT 

C 

IF ( IMIN. GT. ITHRS8)  GO  TO  170 
INTBIT  »  0 
C 

C»M*M»MM«**tMMM*M«M*«***M*MM*MM**M*M**MM 

c*  * 

C«  INTEGER  DISPLACEMENT  ASSUMED  AT  THIS  POINT  * 

C*  * 

£**•••**•*••••«••*•*•**••*••••••*••*•*»•***•*••**•******** 

c 

RCOEF(l)  »  0.0 
RCOEF ( 2)  -  0.0 
RCOEF ( 3 )  -  1.0 
RCOEF(4)  -  0.0 
MIN  -  1 
C 

C  TEST  FOR  NO  DISPLACEMENT 
C 

IFdXOFF. EQ.10  .AND.  IYOFF. EO.  10)  GO  TO  210 
CALL  FILL( IBLK, JBLK, I1M1 , J1M1 ,8 , MTAB, 24 ,16 , ' II ' ) 

NBITS  *  tens  +  20 
GO  TO  210 
170  CONTINUE 
C 

^•••••••••••••••••••••••••••••••••••••••••••••••••••••••* 


2*3 


PERFORM  TEST  TO  DETERMINE  IF  SOB-BLOCK  PROCESSING 
MAT  HELP.  TEST  VARIANCE  OF  ERROR. 


IF ( VA. LT. STDERR  .AND.  ABS(AV) .LT.AVGERR)  GO  TO  240 


****  * 


***  ***• 


•  •  •  *  • 
***•  •••  *••• 


•  •  •  •  •  *  * 
«•••  ••••  *•••  * 


RESET  DISPLACEMENT  MAP 

CALL  FELL(IBLK.JBLK,0, 0.16, MTAB, 24,16, '  ’) 

CORRECT  DATA  RATE  FOR  SOB-BLOCK  PROCESSING 
IF ( INTBIT. HJ.l)  NBITS  -  NBITS  -  18  -  4*NPBITS 

IF( INTBIT. BQ.O .AND. (IIOFF.NE.IO  .OR.  ITOFF. NE. 10) ) NBITS- W ITS-  18 

DO  230  11-1.2 

IBLK1  -  ( 11-1) *8 

DO  230  Jl-1.2 

JBLK1  -  ( Jl-1) *8 

I1M1  -  II  -  1 

J1M1  -  J1  -  1 

KC0N8<1)  -  1 

KCON8( 2)  -  1 

KCON8 ( 3 )  -  28 

KCON8(4)  -  28 


START  DATA  TRANSFER  FOR  8  BY  8  BLOCKS 


IF (IB LOCK. BO. 1)  GO  TO  130 

SECTION  FOR  NON- BORDER  BLOCK  TRANSFER 

I  -  IBLK1 
DO  120  K-1,8 
I  -  I  +  1 
J  -  JBLK1 
DO  120  L-1,8 
J  -  J  +  1 
C8(K.L)  -  C(I,J) 

120  CONTINOE 
GO  TO  150 
130  CONTINOE 

SECTION  FOR  BORDER  BLOCK  TRANSFER 

ITOTAL  -  0 
I  -  EBLK1 
DO  140  K-1,8 
M  -  ISTY  +  I 
I  -  I  +  1 
J  -  JBLK1 
DO  140  L-1,8 
N  -  IX  +  J 
J  -  J  +  1 

ITOTAL  -  ITOTAL  +  IA(M.N) 


CALL  FILLdBLK.JBLK,0,0,16,MrAB,24,16, 'RR') 

CALL  LOCATE(LOC,F,21,21,MIN) 

C 

C  QUADRANT  NUMBER  NOV  CONTAINED  IN  MIN 

C 

c - 

c 

C  CALCULATE  THE  NEW  SET  OF  PREDICTOR  COEFFICIENTS 

C  FOR  THE  CURRENT  BLOCK. 

C 

C  THE  FOUR  VALUES  ARE  CONTAINED  IN  ARRAY  RCOEF. 

C 

C - 

c 

C  LOCATION  OF  THE  4  COEFFICIENTS 

C 

C  RCOEF ( 1)  RCOEF  (2) 

C  - - 

C  RCOEF ( 3 )  RCOEFC4) 

C 
C 
C 

C  SET  UP  THE  REGRESSION  PROBLEM 

C  USE  EVERY  OTHER  VALUE  FOR  THE  PREDICTION  PROCESS 

C 

K  -  0 

DO  80  J-1,16,2 
M  -  J  +  HOFF 
DO  80  1*1.16 ,2 
N  -  I  +  IYOFF 
K  *  K  +  1 
DO  70  L-1,4 

70  X(K.L)  -  D(N+QOFF(MIN, L.l)  , M+QOFF(MIN, L, 2) ) 

80  OUT(K)  -  C(I.J) 

C 

C  GENERATE  THE  PREDICTION  COEFFICIENTS  FROM  THE  MODEL 

C 

CALL  COEFGNU.  OUT.  RCOEF.  K) 

C 

C  QUANTIZE  PREDICTOR  COEFFICIENTS  TO  NPBITS  BITS 
C 

DO  90  1-1,4 

90  RCOEF ( I )  *  ( AINT(RCOEF(I) *RLEVEL  +  .5) ) /RLEVEL 
IF(K.NE.129)  GO  TO  110 
C 

C  IF  INVERSE  DOESN'T  EXIST,  THEN  SET  COEFFICIENTS 

C 

100  CONTINUE 
C 

C  INTEGER  DISPLACEMENT  ASSUMED  AT  THIS  POINT 

C  SET  UP  THE  PREDICTOR  COEFFICIENT  VECTOR 

C 

INTBIT  -  0 
RCOEF(l)  *0.0 
RC0EF(2)  -  0.0 
RC0EF(3)  -1.0 
RCOEF (4)  *  0.0 
MIN  *  1 
C 

C  CHECK  IF  WITH  ERROR  -  NO  DISPLACEMENT 

C 

IFdXOFF. BQ.10  .AND.  IYOFF. BQ. 10)  GO  TO  110 
CALL  FILLdBLK,  JBLK.0 ,0,16  ,MTAB,24,16 ,  ’  II ' ) 

NBITS  -  NBTTS  +  18 
110  CONTINUE 
C 

c - 

c 

C  CALCULATE  THE  PREDICTED  IMAGE  BLOCK 

C 

c - 

C 

CALL  PRED  ( C ,  D ,  ERROR ,36,16,  HOFF ,  IYOFF ,  VA,  AV ,  QOFF ,  MIN ,  RCOEF ,  INTB  IT , 
#  NgJTS, RETCOD, ER16 .DIST.LVLARA) 


nnnnnnnnnfsnnn  nnnnnn  n  noon  oopooooooo  opooo 


START  THE  A  IMAGE  BLOCK  TRANSFER  INTO  TEE  0  SOB- 
MATRIX,  THAT  IS  THE  PAST  FRAME  MATRIX. 


ISTX  -  1ST!  -  10 
ISPX  -  ISTX  +  35 
ITT  -  ISTY  -  10 
IPY  -  ITT  +  35 
ISX  -  ISTX 

mx  -  ispx 

1ST  -  ITT 
IHT  «  IPY 

IF(IBLOCI.EQ.O)  GO  TO  50 
IF(ISX.LT.l)  ISX  *  1 
IFQHX.GT.  JSIZE)  I  HI  =  JSIZE 
IF(IST.LT.l)  1ST  -  1 
IF(IHT.GT.ISIZE)  IHT  *  ISIZE 
M  -  0 
N  -  0 

IF(ISTX.LT.l)  M  -  1  -  ISTX 
IF(ITT.LT.l)  N  -  1  -  ITT 
KCON(l)  -  N  +  1 

KCON( 2)  -  M  +  1 

I  »  M 

DO  60  J-ISX, IHX 
I  -  I  +  1 
L  -  N 

DO  60  I* 1ST, IHT 
L  -  L  +  1 
D(L.K)  -  IA(I,J) 

XCON(3)  -  L 

KCON(4)  -  K 


SOLVE  FOR  THE  METRIC  MATRIX-AND  INTEGER  DISP. 
THE  ACTUAL  METRIC  VALDES  ARE  CONTAINED  IN  THE 
F  MATRIX. 


IF ( IB LOCK. EQ . 1 )  CALL  INITI(F,21,21 .0) 

IF( JPLT.BQ.DCALL  INITI(F,21,21 ,256) 

CALL  METRIC( C ,D,F. 36,21, 16, IMIN , LOC , IB LOCK, KCON, ITHRSH, 11,11) 

TEST  IF  ONLY  INTEGER  DISPLACEMENT  -  IF  SO 
GO  TO  IMAGE  PREDICTION  SECTION. 

IXOFF  -  LOCC2)  -  1 
ITOFF  -  LOC(l)  -  1 

CALL  MPL0T(F,XE,21,IPLT, IMAGE, IBLK.JBLK, LOC, NGRD) 

IF( IMIN. LE. ITHRSH)  GO  TO  100 
INTBIT  -  0 

NON-INTEGER  DISPLACEMENT  ASSUMED  AT  THIS  POINT 


NBITS  -  NBITS  +  18  +  4*NPBITS 
INTBIT  -  1 

DETERMINE  IF  PLOTS  REQUIRED 

IF( JPLT.BQ.l  .AND.  IBLOCX.BQ.O  .AND.  IMAGE. GT. 5) 
CALL  MPL0T(F,2E, 21,  IPLT,  IMAGE,  IBLK.JBLK, LOC,  NGRD) 


SECTION  FOR  NON- INTEGER  PORTION  OF  THE  DISPLACE¬ 
MENT  ESTIMATION. 
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DO  260  IB  LX  -  1 , ISTOP 
ISTY  -  ( IBIX-1) *16  +  1 
ISPY  -  ISTY  +  15 
DO  260  JBIX  -  1 , JSTOP 


COPY  THE  REQUIRED  BLOCK  FROM  IMAGE  B  INTO 
MATRIX  C  THE  PRESENT  FRAME  MATRIX. 


1ST!  IS  THE  X  STARTING  POINTER  FOR  A  AND  B 
ISPX  IS  THE  X  ENDING  POINTER  FOR  A  AND  B 
ISTY  IS  THE  Y  STARTING  POINTER  FOR  A  AND  B 
ISPY  IS  THE  Y  ENDING  POINTER  FOR  A  AND  B 


ISTX  -  ( JBLK-1) *16  +  1 

IX  -  ISTX 

ISPX  -  ISTX  +  15 

RESET  INTBIT 

INTBIT:  0  -  INTEGER  ONLY  DISPLACEMENT 
1  -  NON- INTEGER  DISPLACEMENT 

INTBIT  -  0 

TEST  IF  BLOCK  IS  ON  THE  BORDER.  IF  SO  SET  D  MATRIX  TO  THE 
AVERAGE  OF  THE  B  SOB MATRIX.  LB LOCK  IS  A  BORDER  FLAG 
IF  IB  LOCK  -  0  THEN  NON- BORDER  BLOCK 
IF  IB  LOCK  -  1  THEN  BORDER  BLOCK 

IF (IB LX  .BQ.  1  .OR.  IBLK  .BQ.  ISTOP)  GO  TO  20 
IF (JBIX  .BQ.  1  .OR.  JBLK  .BQ.  JSTOP)  GO  TO  20 


SECTION  FOR  NON-BORDER  BLOCKS 


IB LOCK  -  0 
K  -  0 

DO  10  J- ISTX. ISPX 
K  -  K  +  1 
L  -  0 

DO  10  I- ISTY, ISPY 
L  -  L  +  1 
C(L.K)  »  B(I.J) 

GO  TO  40 


SECTION  FOR  BORDER  BLOCKS 


IB LOCK  -  1 

K  -  0 

ITOTAL  -  0 

DO  30  J-ISTX, ISPX 

K  -  K  +  1 

L  -  0 

DO  30  I- ISTY, ISPY 
L  -  L  +  1 

ITOTAL  -  ITOTAL  +  IA(I,J) 

C(L,K)  -  B(I, J) 

IAVG  -  IFIX(FLOAT( ITOTAL) /256 .  +  .5) 
CALL  INITI(D,36 ,36 , IAVG) 

CONTINUE 
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C  NIB ITS  -  RUNNING  TOTAL  FOR  BITS  USED 

C  NIMAGE  -  NUMBER  OF  IMAGES  TO  BE  PROCESSED 

C  SNRSET  -  MINIMUM  OUTPUT  SNS 

C 

NPB ITS  -  10 
VAREST  =  1 .00 
STDER&  *3.0 
AVGERR  =  3.0 
ISIZE  =  192 
JSIZE  *  128 
IPRT  *  1 

IPLT  *  0 

JPLT  -  0 

IPLT  *  0 

NIB  ITS  *  0 
NIMAGE  -  40 
SNRSET  -  24. 

C 

C  END  OF  PROGRAM  CONSTANT  INITIALIZATION 
C 

RLEVEL  -  FLOAT ( 2*  *NPB ITS ) 

NPB  ITS  -  NPB  ITS  +  1 
ITHRS8  -  VAREST* 12 8 
ITHRSH  -  INT( VAREST* 512 .  +  .5) 

ISTOP  -  ISEZE/ 16 
JSTOP  *  JSIZE/ 16 
SIZE  *  FLOAT(ISIZE*JSZZE) 

SSST  *  0.0 
SSET  *  0.0 
SSPT  *  0.0 
TART  -  0.0 
DRTT  *0.0 

CALL  READUA.  ISIZE,  JSIZE.21) 

C 

C  ALLOW  FOR  MORE  THAN  1  ERROR 

C 

C  CALL  ERRSET(208. 1000. 1.0,0) 

C  CALL  ERRSET(253. 1000. 1.0,0) 

C 

C  INITIALIZE  HISTOGRAM  VECTOR 

C 

CALL  INITR(HIST,1 ,238,0.0) 

CALL  INITKLVLARA.  1,128.0) 

CALL  QUANTI(DIST) 

IF ( IPLT. BQ .1  .AND.  IPLT. BQ.l) CALL  PLOTS ( IBUF, 1 ,15) 

C*********************************** **•*•*•****•**•***•*• 

c* 

C*  START  ACTUAL  IMAGE  SEQUENCE  LOOP 

C* 

C******************************************************** 


IF ( IPRT. GE. 1 ) WRITE ( 6 , 1020) NPB ITS , VAREST. STDERR, AVGERR. SNRSET 
DO  500  IMAGE  -  1, NIMAGE 
CALL  READ(B, ISIZE, JSIZE.21) 

IF(IPRT.  GE. 1 ) WRITE( 6 , 1010) FF, IMAGE 
IF( IPRT. GE.2 ) WRITE( 6 , 103 0) 

IF ( IPRT. GE.2) WRITE(6,1040) 


NBITS  -  A  COUNTER  FOR  THE  DATA  RATE 
NBITS  -  0 


INITIALIZE  MOTION  TAB  ARRAY 
CALL  INITI2(MTAB,24,16 ,16448) 


SET  UP  THE  BLOCK  COUNTERS  TO  STEP  THROUGH  THE 
IMAGE  IN  16  BY  16  BLOCKS.  IB  LX  IS  THE  BLOCK 
ROW  NUWER  AND  JBLK  IS  THE  BLOCK  COLUMN  NUMBER 


PCEC 


PROGRAM  PCEC 


C* 

c* 

c* 

c* 

c* 

PROGRAM  BY:  CARL  BOWLING 

LATEST  UPDATE:  JULY  18,1983 

0 

0 

PROGRAM  FOR  2-D  DELTA  FUNCTION  LOCATION  DETERMIN- 

0 

c* 

ATION 

TO  BE  USED  FOR  MOTION  COMPENSATED  IMAGE* 

c* 

c* 

c* 

C* 

c* 

CODING. 

SUBROUTINES  NEEDED: 

0 

READ 

— 

READS  IN  IMAGES 

c* 

INITR 

- 

INITIALIZES  A  REAL  ARRAY  TO  A  CONSTANT 

c* 

INITI 

— 

INITIALIZES  AN  INTEGER  ARRAY  TO  A  CONST* 

c* 

IN  IT  12 

— 

INITIALIZES  AN  INTEGER *2  ARRAY  TO  A 

0 

c* 

CONSTANT 

* 

c* 

LOCATE 

— 

LOCATES  THE  MINIMUM  QUADRANT  VALUE 

• 

c* 

METRIC 

— 

COMPUTES  THE  DIFFERENCE  METRIC 

• 

C* 

PRED 

— 

CALCULATES  THE  PREDICTED  VALUES 

• 

c* 

COEFGN 

— 

PREDICTION  COEFFICIENT  GENERATOR 

* 

c* 

COMTAL 

— 

OUTPUTS  AN  ARRAY  TO  TAPE  IN  COMTAL  FORM* 

C* 

PLOT 

— 

ALL  THE  CAL COMP  DRIVER  SOFTWARE  AND 

0 

c* 

THE  TRANSLATOR  SOFTWARE  FOR  THE  TEK 

0 

c* 

4662  PLOTTER. 

0 

c* 

INVERT 

— 

SCALES  DATA  0-1  IN  INVERTED  ORDER 

0 

c* 

PUR 

— 

DRIVER  FOR  PURJOY  3D  PLOT  ROUTINE 

0 

c* 

c* 

FILL 

• 

KEEPS  TRACK  OF  BLOCKS  WITH  MOTION 

0 

0 

C « 
C 


C*« 

C* 

C* 

C* 

C •« 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

C 

c 

c 

c 

c 

c 

c 

c 

c 


INTEGER *2  IA(  192 , 128) ,B( 192,128) , EST( 192 ,128) ,MTAB(24,16) ,FF 
INTEGERS  ERRSIG(  192 ,128) 

INTEGER  QOFF(4,4 ,2) ,ICON(4) ,C8(8,8) ,08(28.28) ,E8(8,8) ,KC0N8(4) 
INTEGER  P(21.21) ,C(16,16) ,0(36.36) , ERROR ( 16 . 16 ) ,LVLARA(128) 
INTEGER  ER8( 8,8) ,ER16( 16,16) ,L0C(4) .RETCOD 
REAL  1(64.4) ,0tJT(64) ,HIST(258) ,XHIST(258) ,RC0EP(4) ,1E(21,21) 

REAL  11(2,6)  ,NGRD(4)  .SNRARY (50)  , SNRPAR(SO)  , ORTARY ( 50 )  , TARART ( 50 ) 
REAL  DISTU28) 

COMMON/ THRESH/  SIDERR, AVGERR, SNRSET 
COMMON/ RCODE/  RETCOD 

DATA  QOFF/2*0,2*-1 ,2*1,4*0,2*-1 ,2*1, 2 *0,-1 ,2*0,2*-l ,2*0, -1 ,0,2*1, 
A  2*0. 2*1,0/ 

DATA  FF/3084/ 


PROGRAM  CONSTANT  INITIALIZATION  SECTION 


PROGRAM  CONSTANT  KEY 

NFBITS  -  NUMBER  OF  BITS  TO  QUANTIZE  PREDICTOR  COEFFICIENTS 
(ACTUAL  NUMBER  IS  1  GREATER  TO  INCLUDE  SIGN) 

VAREST  -  ESTIMATE  FOR  AVERAGE  ABSOLUTE  DIFFERENCE  BETWEEN 
CONSTANT  FRAMES 

SIDERR  -  ERROR  STANDARD  DEVIATION  THRESHOLD 
AVGERR  -  ERROR  AVERAGE  THRESHOLD 
I SIZE  -  INPUT  HEIGHT  OF  THE  PICTURE  IN  PIXELS 
JSIZE  -  INPUT  WIDTH  OF  THE  PICTURE  IN  PIXELS 
IPRT  -  PRINTER  FLAG 

0:  DON’T  PRINT  DATA 
1:  PRINT  DATA 

2:  PRINT  ALL  DATA  (STEP  BY  STEP) 

IPLT  -  PLOTTER  DIRECTOR  FLAG 

0:  PLOT  WILL  BE  DIRECTED  TO  4662  PLOTTER  VIA  GPIB 
1:  PLOT  WILL  BE  DIRECTED  TO  THE  4025  SCREEN 
NOTE:  MAXE  SURE  THE  REQUIRED  TXTLIBS  ARE  AVAILABLE 
JILT  -  METRIC  PLOT  FLAG 

0:  NO  METRIC  PLOT  DRAWN 
1 :  METRIC  PLOT  WILL  BE  DRAWN 
KPLT  -  PLOTTER  FLAG 

0:  NO  PLOTTING  IS  NEEDED 
1:  SOME  TYPE  OF  PLOTTING  IS  REQUIRED 
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SUBROUTINE  XFORM 


SUBROUTINE  XFORM(C,  V,XI,  H.N.M) 

c«m*«*mmm****m**mmm**m*m*m***m***m**o*< 

c* 

C*  SUBROUTINE  TO  PERFORM  A  LINEAR  TRANSFORMATION 

C* 

C*  TRANSFORMATION  IS  GIVEN  B7: 

C*  (VHIIHH) 

C* 


REAL  C(N,M) ,  V<2,2)  ,H<M,M)  ,T(6,8)  ,XI(N,M) 
C 

C  CALCULATE  (V)*(XI) 

C 

DO  30  L-1,6,2 

11- L 

12- L+1 

DO  10  1-11,12 
NN  -  I-L+l 
DO  10  J-1,8 
T(I.J)«0.0 
DO  10  C-1,2 
13  -  L  +  I  -  1 

T(I,J)  -  T(I, J)  +  V(NN,K)*XI(I3.J) 

10  CONTINUE 
C 

C  CALCULATE  ((V)(XI))*H 

C 

DO  20  I-I1.I2 
DO  20  J-1,8 
C(I,J)-0.0 
DO  20  K-1.8 

C(I.J)-C(I.J)+T(I.I)*H(K,J) 

20  CONTINUE 
30  CONTINUE 
RETURN 


SUBROUTINE  RINTRP 


FUNCTION  RINTRP<Il,I2,X) 

CM*»MM****M*tM«**M*M*MMMM*»M«**M»*MM*M*M 

c* 

C*  FUNCTION  RINTRP 

C* 

C*  PURPOSE: 

C*  THIS  FUNCTION  IS  USE  TO  PERFORM  A  LINEAR 

C*  INTERPOLATION  BETWEEN  INPUT  INTEGER  POINTS 

C* 

C*  EXPLANATION  OF  CALL  VARIABLES 

C*  II  -  LOWER  BOUND  VALUE 

C*  12  -  UPPER  BOUND  VALUE 

C*  X  -  THE  NON- INTEGER  INTERPOLATION 

C* 

(;•••**•*••*•*•••**•**•••••*••••**•*«**••****«•*•*•*•***••, 

C 

II  -  FLOAT(Il) 

Y2  -  FLOAT (12) 

IF(X.BQ.O.O)  RI  -  Y1 
IF(X.EQ.O.O)  RETURN 
RINTRP  «  Y1  +  (Y2  -  Y1)*X 


non  non  non  non 


50  S2  -  S2  +  R<I)*V1<K) 

DO  60  K-l.N 

G(K)  »  V1(K)/(W  +  S2) 

60  RT  ■  RT  +  A(K)*R(K) 

DO  70  K-l.N 
DO  70  L-l.K 

V(K,L)  -  V(K.L)  -  G(K)*V1(L) 

70  V(L,K)  -  V(K, L) 

E  -  S  -  RT 
DO  80  I-l.N 

80  A(K)  *  A(K)  +  G(K)*E 

SHIFT  THE  PAST  VALUE  VECTOR 
DO  90  L-2.N 

90  R(N*2-L)  -  R(N+1-L) 

R(l)  -  S 
100  CONTINUE 

DO  110  I-l.N 
110  CM(K.I)  -  A(K) 

RESET  R-REGISTER  TO  INITIAL  SIGNAL  VALUE 

DO  120  J-l.N 
120  R(J)  -  XI 

DO  150  J-l , NSF 
S  -  XV(J) 

RT  -  0.0 

DETERMINE  RESIDUAL  SIGNAL 
DO  130  I-l.N 

L30  RT  -  RT  +  A(X)*R(I) 

E-S-RT 
XIN(I.J)  -  E 

SHIFT  R-REGISTER 

DO  140  1-2. N 
140  R(N+2-K)  -  R(NH-I) 

R(l)  -  S 
150  CONTINUE 
160  CONTINUE 
RETURN 
END 
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SUBROUTINE  COMP 


SUBROUTINE  COMPUIN,  NF,  NSF.N,  CM) 
C 

C****** ******************************** 


c* 

* 

c* 

SUBROUTINE  -  COMP 

0 

c* 

PURPOSE: 

0 

c* 

THIS  SUBROUTINE  IS  USED  TO  PERFORM  ADAPTIVE* 

c* 

HYBRID  PICTURE  CODING  (AHPC)  ON  AN  INPUT 

c* 

ARRAY. 

c* 

c* 

EXPLANATION  OF  CALL  VARIABLES: 

c* 

XIN 

-  THE  INPUT  DATA  AND  RESIDUAL  OUTPUT  MATRIX 

c* 

NF 

-  NUMBER  OF  LINES  IN  THE  INPUT  MATRIX 

c* 

NSF 

-  NUMBER  OF  SAMPLES  PER  FRAME- SAMPLES/ LINE 

c* 

N 

-  NUMBER  OF  PREDICTOR  COEFFICIENTS  TO  BE  USED* 

c* 

c* 

c* 

CM 

-  MATRIX  CONTAINING  PREDICTOR  COEFFICIENTS 

DEFINITION  OF  VARIABLE  TERMS: 

c* 

V 

-  THE  ERROR  COVARIANCE  MATRIX 

c* 

A 

-  THE  PREDICTOR  COEFFICIENT  VECTOR 

c* 

VARI  -  INITIAL  VALUE  FOR  ERROR  COVARIANCE  MATRIX 

c* 

W 

-  VARIANCE  OFFSET 

c* 

XV 

-  THE  INPUT  LINE  TEMPORARY  VECTOR 

c* 

6 

-  THE  GAIN  VECTOR 

c* 

R 

-  THE  PAST  VALUE  VECTOR 

c* 

c* 

c* 

E 

-  THE  ERROR  OR  RESIDUAL  TERM 

EXTERNAL  ROUTINES  REQUIRED: 

c* 

NONE 

c* 

c 


c 

REAL  A(6) ,  VL(6)  .HN(NF.NSF)  ,R(6)  ,XV(2 56)  ,0(6) ,  V(6,6)  ,CM(N,NF) 

C 

C  SET  UP  THE  CONSTANTS  AND  INITIAL  VALUES 

C 

DATA  V/ 36*0.0/, A/1.0, -.5, -.2,. 3,. 4, -.5/ 

XNSF  -  NSF 
W  =  1.0 
VARI  =  100.0 
C 

C  IN  THE  DO  160  LOOP.  I  IS  THE  LINE  NUMBER  (1  -  NF) 

C 

DO  160  I-l.NF 
DO  10  J-l.N 
DO  10  K-l.N 
V(J.K)  -  0.0 
IF(J.BQ.K)  V(J,K)  -  VARI 
10  CONTINUE 
C 

C  XI  IS  THE  FIRST  VALUE  OF  EACH  LINE,  USED  TO  INITIALIZE  THE  R  VECT 

C 

XI  =»  XIN(  1,1) 

C 

C  SET  UP  THE  INPUT  VECTOR  AND  THE  PAST  VALUE  VECTOR 

C 

DO  20  J=1 , NSF 
20  XV(J)  -  XIN(I, J) 

DO  30  J=1,N 
30  R(J)  =  XI 
C 

C  IDENTIFICATION  LOOP  (IDENTIFY  THE  PREDICTOR  COEFFICIENTS) 

C 

C 

DO  100  J-l.NSF 
S  -  XV(J) 

S2  -  0.0 
RT  -  0.0 
DO  40  K-l.N 
VI (X)  »  0. 

DO  40  L-1,N 

40  V1(E)  =  VI (K)  +  V(K,L)*R(L) 

DO  50  K«1,N 


IF(MTAB(I,K) .NE. 16448) NTAB  -  NTAB  +  1 
325  CONTINUE 

IF(IPRT.GE.1)WRITE(6,1110)<MTAB(I.K),K=1,16) 

330  CONTINUE 

TARARY (IMAGE)  -  NTAB 
TART  »  TART  +  NTAB 
IF ( IPRT. GE . 1 ) WRITE ( 6 , 1 1 00 ) 

DRATE  -  FLOAT (NB ITS) /SIZE 
DR IT  -  DRTT  +  DRATE 
SSST  =  SSST  +  SSS 
SSET  >  SSET  +  SSE 
SSPT  -  SSPT  +  SSP 
DRTARY(  IMAGE)  =  DRATE 
NTS  ITS  -  NTBITS  +  NBITS 
TDRATE  -  FLOAT(NTB ITS ) /(FLOAT( IMAGE) *SIZE) 

IF ( IPRT. GE. 1 ) WRITE ( 6 . 1120 )  hB  ITS , DRATE. NTB  ITS , TDRATE 
500  CONTINUE 

SNR  -  10. *AL0G10( SSST/ SSET) 

SNRP  -  10. *AL0G10( SSST/ SSPT) 

SNRARY ( NIMAGE+1 )  -  SNR 
SNRPAR ( NIMAG  E+l )  -  SNRP 
DRTARY ( NIMAGE+1 )  »  DRTT/ FLOAT ( NIMAG E ) 

TARARY ( NIMAGE+1 )  -  TART/ FLOAT (NI MAGE) 

NN  a  NIMAG  E  +  1 
WRITE ( 9) NN 

WRTTE(9) ( SNRARY (I). 1=1, NN) 

WRITE( 9) ( SNRPAR (I ) , 1=1 , NN) 

WRITE ( 9) ( DRTARY ( I ) , 1=1 , NN) 

WRITE(9)  (TARARY(I) ,  1=1  ,NN) 

WRITE(9) (HIST(I) ,1=1,256) 

WR3TE(9) (LVLARA(I) , 1=1,128) 

IF( IPRT. GE. 1 ) WRITE ( 6 . 1080) SNR, SNRP 
1010  FORMAT ( IX, A1 , 'IMAGE  SEQUENCE  NUMBER  ',12) 

1020  FORMAT* IX. 'BITS/ COEF.  =  ' , 12 , '  VAREST  -  ',F5.2, 

#  '  STDERR  -  ',F5.2, '  AVGERR  =  ',F5.2,'  SNRSET  =  '.F5.2) 
1030  FORMAT( / , '  Y  X  QUAD  XDIS  YDIS  PRE-AVG  PRE-STD 

#VG  AFT-STD' ,/) 

1040  FORMAT*  IX,  80  ('-')) 

1050  FORMAT(ZX,2I5 , 14 ,215 .2F10.3 ,317) 

1060  FORMAT( IX, 5 15 .4F10.3 ,217) 

1070  FORMAT ( IX, 2F10 . 1 ) 

1080  FORMAT (IX, 'SIGNAL  TO  NOISE  RATIO  *,2F10.5,'  DB') 

1090  FORMAT(lX.Al) 

1100  FORMAT ( IX,  '  ,  _  ,  _ ') 

mo  format(ix.  'pnsnx.AiT;  M  ’) 

1120  FORMAT (IX, '#BITS  =',I7,'  .  DATA  RATE  ',F6.4,'  TOTAL  =  ' 

#  '  CUMULATIVE  DATA  RATE  '  ,F6.4) 

STOP 

END 

BLOCK  DATA 
C 

C  BLOCK  DATA  FOR  POINTERS  FOR  SUBROUTINE  METRIC 
C  IN  ORDER  THEY  ARE  XSTART, XSTOP, YSTART, YSTOP 

C  MAX  IS  21  AND  MIN  IS  1 
C 

COMMON  /PNTR/IRANGE 
INTEGER  IRANGE(4) / 8, 14, 8 ,14/ 


CHANGE 


SUBROUTINE  CHANGE 


Qmmmt 

SUBROUTINE  CHANGEU,  B,  NBY,  N.  IDIR) 

|lt 

c* 

• 

c* 

SUBROUTINE:  CHANGE 

* 

c* 

PURPOSE: 

• 

c* 

THIS  SUBROUTINE  IS  USED  TO  CONVERT 

A  REAL*4* 

c* 

INTO  A  REAL *8  AND  VIS- VERSA  ACCORDING  TO 

c* 

THE  VALUE  OF  IDIR 

c* 

c* 

EXPLANATION  OF  CALL  VARIABLES: 

c* 

A  -  B  -  REAL *4  ARRAY 

c* 

B  -  B  -  REAL *8  ARRAY 

c* 

NBY  -  I  -  1ST  SIZE  VARIABLE  FOR 

A  AND  B 

c* 

N  -  I  -  2ND  <?IZE  VARIABLE  FOR 

A  AND  B 

c* 

IDIR  -  I  -  DIRECTION  FLAG  FOR  TRANSFER 

c* 

+1  MOVE  A  ->  B 

c* 

-1  MOVE  B  ->  A 

c* 

c* 

SUBROUTINES  NEEDED: 

c* 

NONE 

c* 

c***« 

REAL  A(JBY.N) 

REAL *8  B(NBY,N) 

IF(IDIR.LT.O)  GO  TO  20 

DO  10  J«1,N 

DO  10  1=1, NBY 

10 

20 


B(I.J) 
RETURN 
DO  30  J* 


A(I.J) 
1  »N 


COEFGN 


SUBROUTINE  COEFGN 

SUBROUTINE  COEFGN(X, OUT, COEF,  IERN) 

C‘‘‘ . . . .  . 

C* 

C‘ 

C‘ 

C* 

C* 

C‘ 

c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 

REAL  1(64,4) ,XT(4.64) ,XTX(4,4) ,1111(4,4) ,COEF(4) ,OUT(64) ,XTOUT(4) 
REAL  A(4) ,XTXIA(4) ,RCOEF(4) ,XTXIXT(4,64) 

REAL *8  XTX8(4.4),XTXI8(4.4),WORK(100) 

DATA  A/ 4* 1 . / 

CALL  TRANSP(X,XT,64,4.64,4) 

CALL  NMUL(XT,X, XTX, 4 ,64,4 ,4 ,64,4) 

CALL  CHANGE(XTX, XTX8.4 ,4,1) 

CALL  LINV2F (XTX 8 ,4,4,  XTXI8 , 3 , YORK, IERN) 

C 

C  CHECK  FOR  PROBLEMS  WITH  THE  INVERSE  EXISTING,  IF  PROBLEM 

C  DOES  EXIST  RETURN  TO  MAIN  AND  CORRECT 
C 

IF( IERN. BQ. 12 9) RETURN 

CALL  CHANGE( XTXI , XTXI8 , 4 , 4 , -1 ) 

CALL  MMUL(XT, OUT, XTOUT.4 ,64,1,4,64,1) 

CALL  MMUL(XTXI, XTOUT, COEF, 4, 4, 1,4, 4,1) 

C 

C  ADD  CORRECTION  FOR  EQUALITY  AND  INEQUALITY  CONSTRAINTS 

C 

CTOTAL  -  0.0 
TOTAL  -0.0 
DO  10  1-1,4 

CTOTAL  -  CTOTAL  +  COEF(I) 

DO  10  J-1,4 

10  TOTAL  -  TOTAL  +  XTXI(I.J) 

TOTAL  -  1. /TOTAL 

CTOTAL  -  1.  -  CTOTAL 

CALL  MMUL(XTXI, A, XTXIA.4 ,4 ,1 ,4 ,4 ,1) 

TOTAL  -  TOTAL ‘CTOTAL 
DO  20  1-1,4 

20  RCOEF(I)  -  XTXI A ( I ) ‘TOTAL  +  COEF(I) 

RETURN 

END 
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SUBROUTINE:  COEFGN 
PURPOSE: 

THIS  SUBROUTINE  IS  USED  TO  COMPUTE  A  SET  OF 
REGRESSION  COEFFICIENTS  TO  FIT  THE  DATA 
INPUT  IN  THE  MATRIX  X  AND  THE  VECTOR  OUT. 
THE  COEFFICIENTS  ARE  THEN  OUTPUT  IN  THE 
RCOEF  VECTOR  OR  THE  COEF  VECTOR. 

EXPLANATION  OF  CALL  VARIABLES: 

X  -  I  -  MATRIX  OF  SIZE  4  BY  64  WHICH 

CONTAINS  THE  REGRESSION  DATA 
VECTOR  OF  SIZE  64  CONTAINING 
THE  REGRESSION  DATA 
OUTPUT  VECTOR  OF  SIZE  4  THAT 
CONTAINS  THE  INEQUALITY 
CORRECTED  REGRESSION 
COEFFICIENTS. 

OUTPUT  VECTOR  OF  SIZE  4  THAT 
CONTAINS  THE  REGULAR 
REGRESSION  COEFFICIENTS. 

ERROR  CHECK  FOR  NON-EXI STAN CE 
OF  INVERSE 

NOTE:  ONLY  1  OF  THE  COEFFICIENT  GENERATION  METHODS 
CAN  BE  USED.  THE  CHOICE  IS  MADE  BY  PLACING 
THE  NAME  RCOEF  OR  COEF  IN  THE  SUBROUTINE 
STATEMENT  -  THIS  VECTOR  IS  THEN  RETURNED  TO 
THE  CALLING  PROGRAM. 

EXPLANATION  OF  CALL  VARIABLES: 

SUBROUTINES  NEEDED: 

TRANSP  -  TRANSPOSES  A  MATRIX 
MMUL  -  MATRIX  MULTIPLICATION  ROUTINE 
CHANGE  -  CHANGES  REAL *8  TO  REAL ‘4  AND  BACK 
LINV2F  -  MATRIX  INVERSE  (IMSL  ROUTINE) 


OUT 

-  I  - 

RCOEF 

-  0  - 

COEF 

-  0  - 

IERN 

-  0  - 

COMTAL 


C*< 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C*« 


SUBROUTINE  COMTAL 
SUBROUTINE  COMTAL!  IN,  NBY, N,  LOOT) 


SUBROUTINE:  COMTAL 
PURPOSE: 

THIS  SUBROUTINE  IS  USED  TO  TRANSLATE  AN 
INTEGER* 2  ARRAY  TO  1  BYTE  PIXEL  VALUES  AND 
WRITE  IT  OUT  IN  COMTAL  FORMAT  TO  DISK  OR 
TAPE 

EXPLANATION  OF  CALL  VARIABLES: 


IN 
NBY 
N 

LOST 


I 

I 

I 

I 


INTEGER *2  IMAGE  ARRAY 
1ST  SIZE  VARIABLE  FOR  IN 
2ND  SIZE  VARIABLE  FOR  IN 
LOGICAL  OUTPUT  NUMBER  FOR  THE 
DATA  TO  BE  WRITTEN 


SUBROUTINES  NEEDED: 
NONE 


INTEGER*2  IN ( NBY, N) , IWORK( 512) 
LOGICAL*l  OUT(2,512) 

EQUIVALENCE  (IWORK(l) ,OUT(l,l) ) 
DO  20  1=1, NBY 
DO  10  J=1,N 
10  IWORK(J)  =  IN(I.J) 

20  WRITE! LOUT, 30) (0UT(2,J) , J=1,N) 
30  FORMAT!  4  (128A1)) 

RETURN 


SUBROUTINE  FILL 


C*« 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C • 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

C* 

c* 
c* 
c* 
c* 
c* 
c •« 


SUBROUTINE  FILL(  IBLK, JBLK, IB LIS , JBLK8 , IBLKSZ, MDATA, NBY, N, KIND) 


SUB  ROUTINE :  FILL 
PURPOSE: 

THIS  SUBROUTINE  IS  USED  TO  KEEP  AN  OUTPUT 
ARRAY  WITH  THE  TYPE  OF  MOTION  INVOLVED. 

AN  'I'  IN  IHE  OUTPUT  ARRAY  MEANS  THAT  THE 
MOTION  FOR  THAT  BLOCK  OR  SUB- BLOCK  VAS 
INTEGER  ONLY.  AN  'R'  IMPLIES  THAN  SOME 
NON- INTEGER  PORTION  OF  THE  DISPLACEMENT  WAS 
USED 


EXPLANATION  OF  CALL  VARIABLES: 


IB  LX 
JBLK 
IBLK8 
JBLK8 
IBLKSZ 
MDATA 
NBY 
N 

KIND 


I 

I 

I 

I 

I 

0 

I 

I 

I 


Y  LOCATION  OF  MAJOR  BLOCK 
X  LOCATION  OF  MAJOR  BLOCK 

Y  LOCATION  0  SUB-BLOCK 
X  LOCATION  OF  SUB- BLOCK 
CURRENT  BLOCKS IZE 
ARRAY  FOR  INFORMATION 

1ST  SIZE  VARIABLE  FOR  MDATA 
2ND  SIZE  VARIABLE  FOR  MDATA 
MOTION  TYPE  IN  QUOTES : 

'RR'  ->  NON- INTEGER  DISP. 

'II'  ->  INTEGER  DISPLACEMENT 
'  '  ->  NO  DISPLACEMENT 


SUBROUTINES  NEEDED: 
NONE 


INTEGER* 2  MDATA(fBY.N)  .KIND 
MX1  -  2*(IBLK-1)  +  1 
MX2  -  MX1  +  1 
MY1  -  2* (JBLK- 1 )  +  1 
MY2  »  MY1  +  1 
IF ( IBLKSZ. EQ. 8)  GO  TO  10 
MDATA (MX1.MY1)  -  KIND 
MDATA (MX1.MY2)  -  KIND 
MDATA(MX2 . MY1 )  -  KIND 
MDATA(MX2.MY2)  -  KIND 
RETURN 
10  CONTINUE 

MX2  ■  MX1  +  IBLK8 
MY2  -  MY1  +  JBLK8 
MDATA( MX2 . MY2 )  -  KIND 
RETURN 


SUBROUTINE  INITI 


SUBROUTINE  INITI  (A,  N1  ,N2 ,  IV  AL) 

. . 


c* 

c* 

SUBROUTINE:  INITI 

c* 

PURPOSE: 

c* 

IBIS  SUBROUTINE  IS  USED  TO  INITIALIZE  A 

c* 

INTEGER  ARRAY  TO  A  CONSTANT  VALUE 

c* 

c* 

EXPLANATION  OP  CALL 

VARIABLES : 

c* 

A  -  0  - 

INTEGER  ARRAY  TO  INITIALIZE 

c* 

N1  -  I  - 

1ST  SIZE  VARIABLE  FOR  A 

c* 

N2  -  I  - 

2ND  SIZE  VARIABLE  FOR  A 

c* 

IV  AL  -  I  - 

VALUE  THE  ARRAY  IS  INITIALIZED* 

c* 

TO 

• 

c* 

SUBROUTINES  NEEDED: 

• 

c* 

NONE 

* 

c* 

* 

INITIGER  A(N1,N2) 
DO  10  J-1.N2 
DO  10  I-l.Nl 
10  A(I.J)  -  IV AL 


IN IT  12 


SUBROUTINE  IN  IT  12 


SUBROUTINE  INITI2  ( A,  N1  ,N2 ,  IV AL) 


c* 

c* 

SUBROUTINE:  INITI2 

c* 

PURPOSE: 

c* 

THIS  SUBROUTINE  IS  USED  TO  INITIALIZE  A 

c* 

INTEGER* 2  ARRAY  TO  A  CONSTANT  VALUE 

c* 

c* 

EXPLANATION  OF  CALL 

VARIABLES: 

c* 

A  -  0  - 

INTEGER  ARRAY  TO  INITIALIZE 

c* 

N1  -  I  - 

1ST  SIZE  VARIABLE  FOR  A 

c* 

N2  -  I  - 

2ND  SIZE  VARIABLE  FOR  A 

c* 

IV  AL  -  I  - 

VALUE  THE  ARRAY  IS  INITIALIZED* 

c* 

TO 

* 

c* 

SUBROUTINES  NEEDED: 

* 

c* 

NONE 

* 

c* 

* 

INTEGER* 2  A(N1,N2) 
DO  10  J«1,N2 
DO  10  I-l.Nl 
10  A(I, J)  -  IV AL 


1NZTR 


SUBROUTINE  INITR 


SUBROUTINE  INITR(A,N1  ,N2.  VAL) 


SUBROUTINE:  INITR 
PURPOSE* 

'this  SUBROUTINE  IS  USED  TO  INITIALIZE  A 
REAL  ARRAY  TO  A  CONSTANT  VALUE 

EXPLANATION  OF  CALL  VARIABLES: 

A  -  0  -  REAL  ARRAY  TO  BE  INITIALIZED 

N1  -  I  -  1ST  SIZE  VARIABLE  FOR  A 

N2  -  I  -  2ND  SIZE  VARIABLE  FOR  A 

VAL  -  I  -  VALUE  THE  ARRAY  IS  INITIALIZED 

TO 

SUBROUTINES  NEEDED: 

NONE 


REAL  A(N1,N2) 
DO  10  J-1.N2 
DO  10  I-l.Nl 
10  A(I.J)  -  VAL 
RETURN 
END 


INVERT 


SUBROUTINE  INVERT 
SUBROUTINE  INVERT( I A, A,  N) 

(:****••*•*•••*••*•*•*•••••*••••**•**•••***•****•*••*«••**•• . 
c* 

C*  SUBROUTINE :  INVERT 

C*  PURPOSE: 

C*  THIS  SUBROUTINE  IS  USED  TO  INVERT  THE  SCALE 

C*  ORDER  AND  SCALE  THE  INPUT  DATA  TO  BETWEEN 

C*  0  AND  1 

C*  EXPLANATION  OF  CALL  VARIABLES: 

C*  IA  -  I  -  INTEGER  ARRAX  OF  INPUT  DATA 

C*  A  -  0  -  OUTPUT  SCALED  REAL  ARRAY 

C*  N  -  I  -  SIZE  OF  THE  INPUT  AND  OUTPUT 

C*  ARRAYS 

C* 

C*  SUBROUTINES  NEEDED: 

C*  NONE 

C* 

C*MMM****«***M***M«*»*MM*M*M«**M*M***M*M*MM* 

DIMENSION  IA(N.N)  ,A(N.N) 

C 

C  FIND  MAX  AND  MIN  VALUES 
C 

MAX  -  -lOOOOOO 
MIN  -  1000000 

DO  10  J-l.N 
DO  10  I-l.N 

IF(IA(I, J) .GT. MAX) MAX  *  IA(I.J) 

IF(IA(I, J) .LT. MIN) MIN  -  IA(I.J) 

10  CONTINUE 
XMAX  -  MAX 

DIF  -  MAX  -  MIN 
WRITE( 6 . 11 ) MAX. MIN. DIF 

11  FORMAT(  IX,  2 II 0 ,  F10 .4 ) 

C 

C  REORDER  AND  SCALE 

C 

DO  20  J-l.N 
DO  20  I-l.N 

A(I. J)  -  (XMAX  -  FLOAT( IA(I, J) ) ) /DIF 
20  CONTINUE 
RETURN 
END 


LOCATE 


SUBROUTINE  LOCATE 

SUBROUTINE  LOCATE (LOC.  E.  NBY,  N,  MIN) 

(;•*••**•*****•***•**•****»•••**••*••*****••*•' 


c* 

•  i 

c* 

SUBROUTINE :  LOCATE 

•  J 

c* 

PURPOSE: 

• 

c* 

THIS 

SUBROUTINE  IS  USED  FOR  LOCATION  THE 

• 

c* 

NON- INTEGER  PORTION  OF  THE  DISPLACEMENT 

* 

c* 

VIA  FINDING  THE  MINIMUM  QUADRANT 

• 

c* 

• 

c* 

EXPLANATION 

OF  CALL  VARIABLES: 

• 

c* 

LOC 

-  I  -  METRIC  MINIMUM  LOCATION  (X.Y) 

•  j 

c* 

E 

-  I  -  INPUT  MATRIX  HOLDING  METRIC 

•  J 

c* 

NBY 

-  I  -  1ST  SIZE  VARIABLE  FOR  E 

• 

c* 

N 

-  I  -  2ND  SIZE  VARIABLE  FOR  E 

•  • 

c* 

MIN 

-  0  -  MINIMUM  QUADRANT  NUMBER 

c* 

• 

c* 

SUBROUTINES 

NEEDED: 

• 

c* 

NONE 

c* 

* 

c 

INTEGER  Q(4) ,E(NBY,N) , LOC(2) 

COMMON  /PNTR/IRANGE(4) 

C  LOCATE  ABSOLUTE  MINIMUM  QUADRANT  BY  SUMMING  CORNER 

C  VALUES  OF  QUADRANTS  AS  GIVEN  BELOV. 

S  I 

C  IV  I  I 

C  - - - 

C  III  II 

c 

c 

C  TEST  IF  INTEGER  DISPLACEMENT  ESTIMATE  AT  BORDER  OF  METRIC 

C  CALCULATIONS.  (ERROR  WILL  RESULT  IN  THAT  METRIC  VALUES  OUTSIDE 

C  DISPLACEMENT  STEPS  WILL  BE  ZERO. ) 

C 

I  -  LOC(l) 

J  »  LOC (2) 

IM1  =»  I  -  1 
IP1  -  I  +  1 
JM1  *  J  -  1 
JP1  *  J  +  1 

IF(I.BQ.IRANGE(1)  .OR.  I.EQ. IRANGE(2) )  GO  TO  20 
IF(J.BQ.IRANGE<3)  .OR.  J.EQ.IRANGE(4) )  GO  TO  40 
C 

C  DISPLACEMENT  ESTIMATE  NOT  ON  BORDER  OF  METRIC 

C 

Q(l)  -  E(IMl.J)  +  E(IMl.JPl)  +  E(I,JP1) 

Q(2)  *  E(I.JPl)  +  E(IPl.JPl)  +  E(IP1,J) 

Q(3)  =-  E(IPl.J)  +  E(IP1,JM1)  +  E(I,JM1) 

Q( 4)  -  E(I,JM1)  +  E(IMI.JMl)  +  E(IMl.J) 

C 

C  FIND  MIN  QUADRANT 

C 

I  MIN  -  Q  ( 1 ) 

MIN  -  1 
DO  10  M-2.4 

IF(Q(M)  .GE.  IMIN)  GO  TO  10 
MIN  -  M 
IMIN  -  Q(M) 

10  CONTINUE 
RETURN 

20  CONTINUE 

IF(J.BQ.IRANGE<3) .OR. J. EQ. IRANGE(4) )  GO  TO  60 
C 

C  SECTION  FOR  ESTIMATE  ON  BORDER  IN  IHE  Y  DIRECTION 
C 

IF ( I . GQ . IRANGE( 2) )  GO  TO  30 
C 

C  ESTIMATE  LOW  IN  Y  DIRECTION  ONLY 

C  TEST  IF  IN  QUADRANT  II  OR  III 

C 
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Q(2)  -  E(JP1,I)  +  E(IP1, JP1) 

Q<3)  »  E(I,JM1)  +  E(IPl.JMl) 

MIN  -  2 

IF(Q(3) .LT.Q(2))MIN-3 
RETURN 

30  CONTINUE 

HIGH  IN  THE  T  DIRECTION  ONLY 
TEST  IF  IN  QUADRANT  I  OR  IV  ONLY 

Q(l)  -  E(IM1,JP1)  +  E(I.JPl) 

0(4)  -  E(IMl.JMl)  +  E(I,JM1) 

MIN  -  1 

IF(Q(4)  .LT.Q(l) )  MIN  -  4 
RETURN 

40  CONTINUE 

SECTION  FOR  ESTIMATE  ON  BORDER  IN  THE  X  DIRECTION 

IF(  J.BQ.IRANGE(4) )  GO  TO  50 

ESTIMATE  LOW  IN  THE  X  DIRECTION  ONLY 
TEST  IF  IN  QUANDRANT  I  OR  II  ONLY 

Q( 1)  -  E(IMl.J)  +  E(IMl.JPl) 

Q(2)  -  E(IP1,J)  +  E(IPl.JPl) 

MIN  -  1 

IF(Q(2) . LT.Q( 1) )  MIN  -  2 
RETURN 
SO  CONTINUE 

HIGH  IN  THE  X  DIRECTION  ONLY 
TEST  IF  IN  QUADRANT  III  OR  IV  ONLY 


Q(3)  -  E(IPl.JMl)  +  E(IPl.J) 
Q(4)  -  E(I-l.J-l)  +  E(I-l.J) 
MIN  -  3 

IF(Q(4)  .LT.Q(3) )  MIN  -  4 
RETURN 

60  CONTINUE 


SECTION  FOR  ESTIMATE  ON  ONE  OF  THE  CORNERS 

THIS  TILL  FIX  THE  QUADRANT  SUCH  THAT  IT  LIES  ON  THE  INTERIOR 

OF  THE  METRIC  MATRIX. 


MIN  -  1 

IF(I.BQ. IRANGE(2)  .AND.  J. BQ. IRANGE(4) )  MIN  -  2 
IF(I.BQ. IRANGE( 1)  .AND.  J. BQ. IRANGE(4) )  MIN  •  3 
IF(I.BQ.  IRANGE(l)  .AND.  J.  BQ.  IRANGE(2) )  MIN  «  4 
RETURN 
END 
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SUBROUTINE  METRIC 

SUBROUTINE  METRIC(A,  B,D,N1,N2,N3  ,  IMIN.LL,  IB , SCON ,  ITHR,  II ,  J1 ) 


C 


c* 

c • 

SUBROUTINE:  METRIC 

c • 

PURPOSE: 

c* 

THIS  SUBROUTINE  IS  USED  TO  CALCULATE  SOME 

c * 

PREDEFINED 

METRIC  BETWEEN  THE  TWO  INPUT 

c* 

MATRICES  A 

AND  B  AND  OUTPUT  THE  VALUES  OF 

c* 
c * 
c * 

THE  METRIC 

IN  MATRIX  D.  (ALL  INTEGER) 

EXPLANATION  OF  CALL  VARIABLES: 

c* 

A  -  I 

-1ST  MATRIX  FOR  METRIC 

c • 

B  -  I 

-  2ND  MATRIX  FOR  METRIC 

c* 

D  -  0 

-  OUTPUT  MATRIX  CONTAINING  THE 

c* 

METRIC  VALUES  -  METHOD  OF 

c* 

THRESHOLD  EXCEEDING  COUNTING 

c* 

N1  -  I 

-  SIZE  VARIABLE  FOR  B 

c* 

N2  -  I 

-  SIZE  VARIABLE  FOR  D 

c* 

N3  -  I 

-  CURRENT  IMAGE  BLOCKS IZE 

c* 

I  MIN  -  I 

-  MINIMUM  VALUE  FOR  SUM  OF  THE 

c* 

ABSOLUTE  VALUE  OF  THE  ERROR 

c* 

AT  THE  MINIMUM  METRIC  LOCATION* 

c* 

LL( 1)  -  0 

-  OUTPUT  1ST  METRIC  LOCATION 

c* 

LL(  2)  -  0 

-  OUTPUT  2ND  METRIC  LOCATION 

c* 

IB  -  I 

-  BORDER  BLOCK  FLAG  (0/1) 

c* 

KCON  -  I 

-  START/ STOP  FLAG  VECTOR 

c* 

ITHR  -  I 

-  DIFFERENCE  THRESHOLD  FOR  1ST 

c* 

TEST  RETURN 

c* 

11  -  I 

-  A  PRIORI  X  INTEGER  ESTIMATE 

c* 

FOR  THE  DISPLACEMENT 

c* 

J1  -  I 

-  A  PRIORI  Y  INTEGER  ESTIMATE 

c* 

c* 

c* 

FOR  THE  DISPLACEMENT 

SUBROUTINES  NEEDED 

c* 

c* 

NONE 

. . . 

INTEGER  LL(4) ,D(N2,N2) ,A(N3 ,N3> ,B(N1,N1) ,KCON(4) ,IPOS(256,2) 
INTEGERS  E<16.16) 

COMMON  /PNTR/  IRANGE(4) 

C  I MIN  -  2000000 

KTHR  =  3 
KTEST  =  N3  **2 
IISIZE  -  ETEST 
SIZE  =«  IISIZE 
C 

C  SET  UP  INITIAL  DISPLACEMENT  ESTIMATE  (A  PRIORI  GUESS) 

C 

ISTART  -  II 
ISTOP  -  II 
JSTART  -  J1 
JSTOP  -  J1 
I COUNT  -  1 
C 

C  DETERMINE  IF  CURRENT  BLOCK  IS  BORDER  BLOCK 

C 

IF(IB.BQ.l)  GO  TO  80 
C 

(;•••••••••••••••••••••»***••****••••*•«**•**•**•*••••*•••* 

C*  • 

C*  SECTION  FOR  NON-BORDER  BLOCKS  * 

C*  * 

. . . 

c 

10  CONTINUE 

ETEST  -  IISIZE 
C 

C  THE  DO  30  LOOPS  ADJUST  THE  DISPLACEMENT  ESTIMATE 
C 

DO  30  I- ISTART, ISTOP 
IM1  »  I  -  1 
DO  30  J« JSTART. JSTOP 
C 
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KCOUNT  IS  USED  FOR  THE  THRESHOLD  EXCEEDING  COUNTING  METHOD 
ISUM  IS  USED  FOR  SUM  OF  ABSOLUTE  VALUE  OF  DIFFERENCE  METHOD 

ICOUNT  »  0 
KSUM  =*  0 
JM1  -  J  -  1 

THE  DO  20  LOOPS  EVALUATE  A  SINGLE  METRIC  VALUE 

DO  20  K»1  ,N3 
IK  *  K  +  IM1 
DO  20  L  -  1.N3 
JL  *  L  +  JM1 

KDIF  -  IABS(A(L,K)  -  B(JL, IK)) 

KSUM  -  KSUM  -i-  KDIF 
IF ( KDIF . GT.KTHR) KCOUNT  -  K COUNT  +  1 
20  CONTINUE 

D(J.I)  -  KCOUNT 
E(J.I)  »  KSUM 

IF (KCOUNT. GT.KTEST)  GO  TO  30 

RESET  COUNTING  METHOD  POINTERS  AND  MINIMUM  VALUE 

ETEST  *  KCOUNT 
LL(  1)  -  J 
LL(  2)  =  I 
30  CONTINUE 

TEST  FOR  POSSIBLE  MULTIPLE  MINIMA  (NCOUNT  -  #MINIMUMS) 

N COUNT  -  0 

DO  40  I- I START, ISTOP 
DO  40  J-JSTART. JSTOP 
IF(D(I, J) .NE.KTEST)GO  TO  40 
NCOUNT  -  NCOUNT  +  1 
IPOS ( NCOUNT, 1)  -  I 
IPOS (NCOUNT, 2)  -  J 
40  CONTINUE 

IF ( NCOUNT. BQ.l)  GO  TO  70 

MULTIPLE  MINIMA  POINTS  FOUND,  RETRY  WITH  SMALLER  ALLOWABLE  ERROR 

KTHR  =•  KTHR  -  1 
IF(KTHR.GE.l)  GO  TO  10 
50  CONTINUE 

MULTIPLE  MINIMA  POINTS  AT  SMALLEST  ALLOWABLE  ERROR 
TAKE  THE  SMALLEST  DISPLACEMENT  TO  BE  THE  ESTIMATE 

MM IN  -  242 

DO  60  INUM  -  1, NCOUNT 

MOT  -  (IPOS( INUM, 1) -11) *(IP0S( INUM, 1)-11) +(IPOS( INUM, 2) -11) * 

»  (IPOS  (INUM,  2) -11) 

IF(MOT.GT.MMIN)  GO  TO  60 

SET  THE  COUNTERS  TO  THE  SMALLEST  DISPLACEMENT 

MM IN  -  MOT 

LL( 1)  »  IPOS ( INUM, 1) 

LL(2)  *  IP0S( INUM, 2) 

60  CONTINUE 
70  CONTINUE 

IMIN  -  E(LL( 1) ,LL( 2) ) 

IF(IMIN.LE. TIER  .OR.  ICOUNT. BQ. 2)  RETURN 
ICOUNT  -  2 

SET  FULL  RANGE  METRIC  POINTERS 


80 


ISTART  - 
ISTOP  - 
J START  - 
JSTOP  - 
GO  TO  10 
CONTINUE 
ETEST  -  N3*«2 


DLANGE(l) 

OLANGE(2) 

HANGE(3) 

IRANGE(4) 
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C  IMIN  »  2000000 

C 

(;•*•*•*•**•*••*•*•**••**•*••******«*•*•***•**•************ 

C*  * 

C«  SECTION  FOR  BORDER  BLOCKS  * 

C*  * 

. . . . ******************** 

c 

DO  110  I* ISTART, ISTOP 
IM1  -  I  -  1 
DO  110  J*  JSTART, JSTOP 
C 

C  KCOUNT  IS  USED  FOR  IRE  THRESHOLD  EXCEEDING  COUNTING  METHOD 

C  KSUM  IS  USED  FOR  SUM  OF  SQUARES  OF  DIFFERENCE  METHOD 

C  KK  IS  A  COUNTER  FOR  NUJBER  OF  PIXEL  IN  BORDER  BLOCKS 

C 

KCOUNT  -  0 
KSUM  -  0 
KK  -  0 
JM1  =*  J  -  1 
DO  90  K-1.N3 
IK  -  K  +  IM1 
C 

C  DETERMINE  IF  OUTSIDE  OF  BORDER 

C 

IF( IK. LT. KCON( 2)  .OR.  IK.GT.KC0N(4) )  GO  TO  90 
DO  90  L-1.N3 
JL  -  L  +  JM1 
C 

C  DETERMINE  IF  OUTSIDE  OF  BORDER 

C 

IF ( JL. LT. KCON ( 1)  .OR.  JL.GT.KCON( 3) )  GO  TO  90 

KK  -  KK  +  1 

KDIF  -  IABS(A(L, K)  -  B(JL.IK)) 

KSUM  -  KSUM  +  KDIF 
IF(KDIF.GT.KTHR) KCOUNT  -  KCOUNT  +  1 
90  CONTINUE 
XKK  -  KK 

IF(KK.NE.O)  GO  TO  100 
C 

C  NO  PIXEL  OVERLAP 

C 

KSUM  -  255  •  IISIZE 
KCOUNT  -  IISIZE 
KK  -  IISIZE 
100  CONTINUE 
C 

C  AT  LEAST  SOME  PIXEL  OVERLAP 

C 

IF(ICOUNT.BQ.l)  GO  TO  105 
IF(KK.EQ. IISIZE)  GO  TO  105 
C 

C  ADJUST  FOR  NON-FULL  BLOCK 

C 

KSUM  -  IFIX( ( SIZE/FLOAT(KX) ) *KSUM) 

COUNT  -  KCOUNT 

KCOUNT  *  INT ( COUNT* (SIZE/XKK)  +  .5) 

105  D(J.I)  -  KCOUNT 
E( J, I)  *  KSUM 

IF(KCOUNT.GE.KTEST)  GO  TO  110 
C 

C  RESET  COUNTING  METHOD  POINTERS  AND  MINIMUM 

C 

KTEST  -  KCOUNT 
LL(  1)  -  J 
LL(  2)  -  I 

C  IMIN  -  KSUM 

110  CONTINUE 
C 

C  TEST  FOR  POSSIBLE  MULTIPLE  MINIMA  (NCOUNT  -  #MINIMUMS) 

C 

NCOUNT  «  0 

DO  120  I- ISTART. ISTOP 
DO  120  J-JSTART, JSTOP 
IF(D(I, J)  .NE.  KTEST)  GO  TO  120 


nnn  non  nnnn  non 
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NCOUNT  =■  NCOUNT  +  1 
IPOS ( NCOONT, 1 )  =  I 
IPOS ( NCOUNT, 2 )  -  J 
120  CONTINUE 

IF ( NCOUNT. BQ.l)  GO  TO  140 

MULTIPLE  MINIMA  POINTS  FOUND,  RETRY  WITH  SMALLER  ALLOWABLE  ERROR 

KTHR  *  KTHR  -  1 

IF (KTHR  .GE.  1)  GO  TO  80 

MULTIPLE  MINIMA  POINTS  AT  SMALLEST  ALLOWABLE  ERROR 
TAKE  SMALLEST  DISPLACEMENT  TO  BE  DISPLACEMENT 

MM IN  =  242 

DO  130  INUM  -  1, NCOUNT 

MOT  =  ( IPOS ( INUM,  1 )  -  ll)*(IPOS(INUM,  1)  -  11) 

#  ( IPOS ( INUM, 2 )  -  ll)*(IPOS(INUM,2)  -  11) 

IF(MOT.GT.MMIN)  GO  TO  130 

SET  COUNT  TO  SMALLEST  DISPLACEMENT 

MM IN  -  MOT 

LL( 1)  -  IPOS ( INUM, 1) 

LL(2)  =»  IPOS  (INUM,  2) 

130  CONTINUE 
140  CONTINUE 

IMIN  -  E(LL( 1) ,LL(  2) ) 

IF  ( IMIN .  LE.  ITHR  .OR.  I COUNT.  BQ. 2)  RETURN 
I COUNT  -  2 

SET  FULL  RANGE  METRIC  POINTERS 

ISTART  =  IRANGE(l) 

ISTOP  -  IRANGE( 2) 

JSTART  -  IRANGE(3) 

JSTOP  -  IRANGE(4) 

GO  TO  80 
END 
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SUBROUTINE  MMUL(A, B, C, NA,MA,MB,N1 ,N2,N3) 

c«*m«***m***mm*«m***mm««»**m*»***m***m**mm« 

c* 

C*  SUB  ROUTINE  :MMUL 

C*  PURPOSE: 

C*  THIS  SUBROUTINE  IS  USED  TO  MULTIPLY  TWO 

C*  ARRAYS  TO  FORM  A  THIRD  IN  THE  FORMji 

C*  (A) *(B)  =  (C) 

C* 

C*  EXPLANATION  OF  CALL  VARIABLES: 


c* 

A 

— 

I  -  1ST  REAL  INPUT  MATRIX 

c* 

B 

— 

I  -  2ND  REAL  INPUT  MATRIX 

c* 

C 

— 

0  -  REAL  OUTPUT  MATRIX 

c* 

NA 

— 

I  -  1ST  SIZE  VARIABLE  OF  A  AND  C 

c* 

MA 

— 

I  -  2ND  SIZE  VARIABLE  OF  A,  1ST  B 

c* 

MB 

— 

I  -  2ND  SIZE  VARIABLE  OF  B  AND  C 

c* 

N1 

— 

I  -  USABLE  FORM  OF  NA 

c* 

N2 

— 

I  -  USABLE  FORM  OF  MA 

c* 

N3 

- 

I  -  USABLE  FORM  OF  MB 

C* 

C*  SUBROUTINES  NEEDED: 

C*  NONE 

C* 

(;*•*••****«•••*•*••«**••*••••*•*•***•• 

REAL  A(NA,MA) ,B(MA,MB) ,C(NA,MB) 
DO  20  I-l.Nl 
DO  20  J-1.N3 
T  -  0. 

DO  10  K-1.N2 

10  T  »  T  +  A(I,K)*B(K, J) 

20  C(I.J)  -  T 


! 

i 
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MPLOT 

SUBROUTINE  MPLOT 

SUBROUTINE  MPLOT(F,XE, ISIZE, IPLT, IMAGE,  IBLK, JBLK, LOC, NGRD) 

C****M*MMM»»**MO**M***M**MM*****MMM**M******** 

c*  * 

C*  SUBROUTINE  MPLOT  * 

C*  * 

C*  PURPOSE:  * 

C*  THIS  SUBROUTINE  IS  USED  TO  PLOT  A  3-D  DATA-* 

C*  BASE  IN  TWO  DIMENSIONS.  IN  THIS  CASE  IT  IS* 

C*  USED  TO  PLOT  DIFFERENCE  METRIC  (ACTUALLY  AN* 

C*  INVERTED  SCALED  VERSION) .  * 

C*  * 


c* 

EXPLANATION  OF 

CALL  VARIABLES: 

c* 

F 

-  I  - 

METRIC  ARRAY  TO  BE  PLOTTED 

c* 

XE 

-  T  - 

REAL  WORK  ARRAY  VERSION  OF  F 

c* 

ISIZE 

-  I  - 

SIZE  OF  SQUARE  METRIC 

c* 

IPLT 

-  I  - 

DISPLAY  DEVICE  FLAG 

c* 

-  0 

PLOT  TO  4662  PLOTTER 

c* 

-  1 

PLOT  ON  4025  DISPLAY 

c* 

IMAGE 

-  I  - 

INTEGER  VALUE  OF  FRAME  NUMBER 

c* 

IBLK 

-  I  - 

BLOCK  LOCATION  VARIABLE 

c* 

JBLK 

-  I  - 

BLOCK  LOCATION  VARIABLE 

c* 

LOC 

-  I  - 

MINIMUM  VALUE  LOCATION 

c* 

NGRD 

-  I  - 

ARRAY  RBQUIRED  BY  THE  3D  PLOT 

C*  * 

C*  SUBROUTINES  NEEDED:  * 

C*  INVERT  -  CHANGES  ORDER  AND  SCALE  0-1  * 

C*  ERASE  -  CLEARS  4025  SCREEN  * 

C*  FACTOR  -  SCALES  DATA  BEFORE  PLOTTING  * 

C*  PUR  -  3-D  PLOT  ROUTINE  DRIVER  * 

C*  PLOTS  -  CAL  COMP  PLOTTING  SOFTWARE  PACKAGE  • 

C*  • 

C*********************************************************** 

INTEGER  F(  ISIZE,  ISIZE).  LOC  (4) 

REAL  KEdSIZE,  ISIZE)  ,NSRD(4)  ,XY(2,6) 

LUNIT  -  6 
C 

CALL  INVERT(F,XE,  ISIZE) 

C 

C  SET  3D  PLOT  PARAMETERS  -  LOOK  AT  SUBROUTINE  PUR  FOR  MORE  DETAILS 

C 

NGRD(l)  -  0 
NGRD (2)  -  1 
NGRD(3)  «  1 
NGRD (4)  -  1 
X  -  0.0 
1-5.2 

IF ( IPLT.  BQ  .1 )  CALL  ERASE 
IF( IPLT. BQ.O) CALL  PLOT(X, Y,-3) 

IF (IPLT. BQ.O)  CALL  PL0TS(IBUF,1,15) 

IF( IPLT. BO. 1) CALL  FACTOR( .400) 

IF( IPLT. BQ.O) CALL  FACTOR( .40) 

CALL  PUR(XE,  ISIZE,  ISIZE,  ISIZE,XT,0,NGRD,15) 

IF ( IPLT. BQ . 1 ) CALL  FACTOR(.40) 

IF (IPLT. BQ.O) CALL  FACTOR( .40) 

CALL  PL0T(0.0,4 .3,-3) 

C  CALL  SYMBOL( 1. ,6 .9, .20, 'FRAME' ,0.0,5) 

C  CALL  SYMBOL ( 1. ,6 .6 , .20, ' IBLOCK' ,0 . ,6) 

C  CALL  SYMB0L(1. ,6 .3 , .20, ' JBLOCK' ,0. ,6) 

XFRAME  -  IMAGE  +  1 
XBLOCK  -  JBLK 
YBLOCK  -  IBLK 

C  CALL  NUM)ER( 2.6 ,6 .9 , .20 ,XFRAME, 0 .0 ,-l) 

C  CALL  NUM3ER(2.6 ,6 .6 , .20, XBLOCK, 0 .0,-1) 

C  CALL  NU1BER(2.6 ,6 .3 , .20 , YBLOCK, 0 .0 ,-l) 

C  CALL  SYMBOLd.  ,6 . ,  .20, 'X  SHIFT', 0., 7) 

C  CALL  SYMBOL(l. ,5 .70, .20, 'Y  SHIFT',0.,7) 

XSHIFT  -  LOC(l)  -  11 
YSHIFT  -  LOC (2)  -  11 

C  CALL  NUMBER(2. 6, 6.0, .20, XSHIFT, 0. 0,-1) 

C  CALL  NUMBER(2 .6 ,5.7, .20 , YSHIFT, 0 .0,-1) 

C  IF(IPLT.BQ.l)  CALL  SYMBOL( 0.0 ,0 .0 . .10, ' . ' ,0 .0 ,1) 

IF(IPLT.BQ.l)  CALL  AN MODE 
IF ( IPLT. BQ . 1)  CALL  TSEND 
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IF(IPLT. BQ.O)  CALL  PLOT<10. ,10. ,999) 

WRITE (LUNIT, 400) 

BP. An  DC MBIT  ARGUMENT  FOR  PLOTTER  DELAY 
READ (5. 410) INSWER 

400  FORMAT (IX.  '/‘INPUT  ANY  SINGLE  DIGIT  NUMBER  TO  CONTINUE') 
410  FORMAT(Il) 

RETURN 

END 


SUBROUTINE  PEED 


C*« 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

C*« 


SUBROUTINE  PRED(C,D,  ES.  I  CD,  IS,  HF,  ITF,  STD,  AVG,  QOFF,  MIN.  COEF,  IFLAG, 
#  NBITS , SETCOD, F, DIST, LVLARA) 


SUB  ROUTINE:  PEED 
PURPOSE: 

THIS  SUBROUTINE  IS  USED  TO  CALCULATE  THE 
CURRENT  IMAGE  FRAME  FROM  INFORMATION 
CONTAINED  IN  THE  PREVIOUS  IMAGE  FRAME  AND 
THE  PREDICTION  COEFICIENTS  GENERATED  FOR 
THE  CURRENT  FRAME 


EXPLANATION  OF 
C 


CALL 
-  I  - 


-  I  - 


ES 


-  0  - 


I  CD 

IS 

IXF 

ITF 

STD 

AVG 

QOFF 

MIN 

COEF 

IFLAG 

NBITS 

RETCOD 

F 

DIST 


I 

I 

I 

I 

0 

0 

I 

I 

I 

I 

0 

0 

0 

I 


VARIABLES: 

INTEGER  AREA?  FOR  CURRENT 
IMAGE  BLOCK 

INTEGER  ARRAX  FOR  PREVIOUS 

ESTIMATED  IMAGE 

INTEGER  ARRAY  FOR  CURRENT 

ESTIMATED  IMAGE 

SIZE  OF  THE  D  MATRIX 

SIZE  OF  THE  ES  AND  C  MATRICES 

X  OFFSET  VARIABLE 

Y  OFFSET  VARIABLE 

STANDARD  DEV  OF  THE  ERROR 

AVERAGE  OF  THE  ERROR 

ARRAY  OF  POINTER  VALUES 

MINIMJM  QUANDRANT  NUJBER 

PREDICTION  COEFFICIENT  VECTOR 

INTEGER  DISPLACEMENT  FLAG 

DATA  RATE  COUNTER 

LARGE  ERROR/ BLKSIZE  BIT 

OUTPUT  YORK  ARRAY  (ERROR?) 

PREDICTOR  GAIN  (IN  DB) 


SUBROUTINES  NEEDED: 

QUANTZ  -  SETS  VALUES  FOR  ADAPTIVE  QUANTIZER 


INTEGER  C(IS, IS) ,D(ICD, ICD) ,Q0FF(4,4,2) ,ES(IS. IS) , RETCOD 
INTEGER  F(IS, IS) ,G(lti,16) , LVLARA(l) 

REAL  XLEVEL(128) ,XOUT(128) ,COEF(4) ,DIST(1) 

COMMON/ THRESH/  STDERR, AVGERR, SNRSET 

SIZE  -  IS**2 

RETCOD  -  0 

KSTART  -  1 

IS  TOP  -  4 


CHECK  FOR  INTEGER  DISPLACEMENT 


10 


IF ( IFLAG. NE.O)  GO  TO  10 

KSTART  -  3 

IS  TOP  -  3 

KA  -  0 

NA  -  0 

KV  -  0 

NV  -  0 

DO  30  J-l.IS 
JJ  -  J  +  IXF 
DO  30  I-l.IS 
II  -  I  +  IYF 
TEMP  “  0.0 


MAKE  PREDICTION  BASED  ON  PREDICTION  COEFFICIENTS 


20 


DO  20  RESTART,  KSTOP 

TEMP  -  TEMP  +  COEF(M)*FLOAX(D(II+QOFF(MIN,M,l) , JJ+QOFF(MIN, M,2) ) ) 
CONTINUE 

ES(I,J)  -  TEMP  +  .5 
TERROR  -  C(I,J)  -  ES(I,J) 

F(I, J)  -  IERROR 
KA  -  KA  +  IERROR 
KV  -  KV  +  IERROR* IERROR 
NV  -  NV  +  C(I, J)*C(I,J) 


onn  nnnnn  n  nnnnnn  nnnnn  non  nnnn  nnnnn 


PRED 


NA  *  NA  +  C(I,J) 

30  CONTINUE 
AVGK  -  XA 
AV6N  -  NA 
VARX  -  KV 
VARN  -  NV 

VARXX  »  (VARX  -  (AVGK*AVGK/SIZE) )/ (SIZE-1 . ) 

VARNN  *  ( VARN  -  ( AVGN*AVGN/ SIZE) ) / ( SIZE-1 . ) 

IF (VARNN. EQ.O .0) SNR  =  SNRSET 
IF( VARNN. BQ. 0.0) GO  TO  15 
SNR  -  10. *AL0G10( VARNN/ VARKX) 

15  STD  -  SORT (VARXX) 

AVG  -  AVGK/ SIZE 

IF(STD.LT.  STDERR  .AND.  ABS(AVG)  .LT.AVGERR)  GO  TO  70 

RETURN  AND  RETRY  IF: 

1)  LARGE  ERROR  FOR  16X16  BLOCK 

2)  LARGE  ERROR  FOR  INTEGER  DISPLACEMENT  ON  8X8  BLOCK 


IFUS.BQ.16)  RETURN 
IF  ( IFLAG .  BQ .  0 )  RETCOD  =  1 
IF (IFL AD. EQ.O) RETURN 

CORRECT  FOR  LARGE  ERROR  WITH  BLOCK  ADAPTIVE  QUANTIZER 
DETERMINE  QUANTIZER  GAIN  REQUIRED  TO  OBTAIN  SET  SNR 

RBQGAN  -  SNRSET  -  SNR 

QUANTIZE  THE  MEAN  AND  VARIANCE  BEFORE  ERROR  QUANTIZATION 

AVG1  -  (FLOAT( INT( AVG*2 .0  +  SIGN( .5 , AVG) ) ) ) /2.0 
STD1  -  ( FLOAT( INT( STD*8 . 0  +  .5)))/8.0 

CHECK  FOR  QUANTIZER  LEVEL  OVERFLOW 

IF(ABS(AVG1) .GT.32.0) VRITE(6.55) 

IF(  STD1 .  GT.  64 . )  WRITE(  6 . 56  ) 

55  F0RMAT(1X, 'AVERAGE  ERROR  TO  LARGE.  LOOK  AT  PRED. ') 

56  FORMAT ( IX. ’ STANDARD  DEVIATION  ERROR  TO  LARGE,  LOOK  AT  PRED.') 
IF(ABS( AVG1 ) .GT.3  2 .0 ) AVG1-32 . *SIGN( 1 . , AVG) 

IF( STD1 .GT. 64 .0 ) STD1-64 .0 

SUBROUTINE  QUANTZ  WILL  DETERMINE  THE  NUJfflER  OF  LEVELS  REQUIRED 
RBQGAN  -  10.*AL0G10( VARKK/3 .) 

CALL  QUANTZ  (RBQGAN,  AVG1 ,  STD1 ,  LEVEL,  XL EVEL,  XOUT,  DIST) 

I® ITS  -  I® ITS  +  23  +  INT( SIZE* ALOG (FL0AT( LEVEL) )/AL0G( 2.)  +  1.0) 
LVLARA( LEVEL)  *  LVLARA( LEVEL)  +  1 

IF ( LEVEL . LE . 1 28 ) WRITE ( 6 . 11 ) LEVEL , RBQGAN , AVG1 , STD1 , VARNN, VARKK, AVGK 
11  FORMAT( IX, 'REQUIRES’, 13,'  BITS', 6F10. 4) 

7  BITS  FOR  THE  MEAN 

9  BITS  FOR  THE  STD 

7  BITS  FOR  THE  NUMBER  OF  LEVELS 

DO  60  J-l.IS 
DO  60  I-l.IS 
IER  -  C(I,J)  -  ES(I.J) 

X  »  IEE 

DO  40  K-l, LEVEL 

IF(X. GT.XLEVEL(K) )  GO  TO  40 

IQUAN  -  INT(X0UT(K)  +  SIGN( .5 ,XOUT(K) ) ) 

GO  TO  50 
40  CONTINUE 
50  CONTINUE 

ES(I.J)  -  ES(I.J)  +  IQUAN 
0(1, J)  -  C(I,J>  -  ES(I,J) 

60  CONTINUE 
70  CONTINUE 

RESTRICT  OUTPUT  FOR  0-255. 

DO  80  J-l.IS 
DO  80  I-l.IS 

IF(ES(I, J) .GT.255)  ES(I,J)  -  255 
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SUBROUTINE  TRANSP 
SUBROUTINE  TRANSP  ( A,  B,  Nl,  N2,  N3  ,N4) 

c* 

C*  SUBROUTINE : TRANSP 

C*  PURPOSE: 

C*  THIS  SUBROUTINE  IS  USED  TO  TRANSPOSE  ARRAYS 

C* 

C*  EXPLANATION  OF  CALL  VARIABLES: 


c* 

A 

-  I  - 

INPUT  DATA  ARRAY 

c* 

B 

-  0  - 

OUTPUT  REAL  DATA  ARRAY 

c* 

Nl 

-  I  - 

1ST  ACTUAL  SIZE  VARIABLE 

OF 

A 

c* 

IN  CALLING  PROGRAM 

c* 

N2 

-  I  - 

2ND  ACTUAL  SIZE  VARIABLE 

OF 

A 

c* 

IN  CALLING  PROGRAM 

c* 

N3 

-  I  - 

1ST  USING  SIZE  VARIABLE 

c* 

N4 

-  I  - 

2ND  USING  SIZE  VARIABLE 

c* 

c* 

SUBROUTINES 

NEEDED: 

c* 

NONE 

c* 

Qmm+i 

REAL  A(N1,N2) 

,B(N2 ,N1) 

DO  10  I-1.N3 

DO  10  J-1.N4 
10  B( J, I)  -  A(I.J) 
RETURN 


QUANTI 


SDB ROUTINE  QUANTI 
SUBROUTINE  QUANTI  (DIST) 

(;•**•••***••*•**•*••***••*•*••***•*•«**•*«*••**••••«•*•***** 

C*  * 

C*  SUBROUTINE  QUANTI  * 

C*  PURPOSE:  * 

C*  THIS  SUBROUTINE  IS  USED  TO  GENERATE  THE  * 

C*  FUNCTIONAL  VALUES  FOR  A  NORMAL(O.l)  * 

C*  DISTRIBUTION.  * 

C*  * 

C*  EXPLANATION  OF  CALL  VARIABLES  :  * 

C*  DIST  -  0  -  ARRAY  CONTAINING  GENERATED  * 

C*  DATA  • 

C*  * 

C*  SUBROUTINES  NEEDED:  • 

C*  MDNRIS  -  GENERATE  VALUES  FOR  GOASSIAN  DIST  • 

C*  * 

(:••********•****•*•****•••••••*•••*•*•«•*•**••**••••***•**•* 

REAL  X(lOOl) ,F(1001) ,0UT(1003) ,XOUT(128) ,XLEVEL(128) .DIST(l) 
C 

C  GENERATE  THE  VALUES  FOR  THE  N(O.l)  DISTRIBUTION 
C 

DO  10  1*1,1001 

X(I)  *  (FLOAT(I)  -  S01J/100. 

F(I)  -  ,398942*EXP(-.5*(X(I)*X(I) ) ) 

10  CONTINUE 

DIST(l)  =  0. 

DO  90  LEVEL*2 ,128 
XLEVL  *  LEVEL 
START  -  1./ XLEVL 
NLEVL1  *  LEVEL  -  1 
PERCNT  *  START 
DO  20  1*1 , NLEVL1 

CALL  MDNRIS  (PERCNT,  XLEVEL(  I),  IER) 

20  PERCNT  -  PERCNT  +  START 
XLEVEL( LEVEL)  -  1000. 

BEGIN  *  START/2. 

PERCNT  -  BEGIN 
DO  30  1*1, LEVEL 

CALL  MDNRIS  ( PERCNT,  XOUT(  I),  IER) 

30  PERCNT  -  PERCNT  +  START 
DO  40  1*1,1001 
40  OUT(I)  -  0.0 
DO  70  1*1,1001 
VALUE-  (FLOAT( 1-501 ) ) /100. 

C - DETERMINE  WHICH  LEVEL  IT  FALLS  WITHIN 

DO  50  11=1, LEVEL 
IF( VALUE. GT.XLEVEL( II))  GO  TO  50 
K  -  II 
GO  TO  60 
50  CONTINUE 
60  CONTINUE 

C - XOUT(K)  IS  THE  AMOUNT  OF  SHIFT 

IV AL  -  INT( (VALUE  -  XOUT(K))*100  +  501.5) 

OUT(IVAL)  *  OUT(IVAL)  +  F(I) 

70  CONTINUE 

C - CALCULATE  THE  DISTORTION 

DIS  -  0.0 
DO  80  1-1,1000 
XX  -  X(I)*X(I) 

80  DIS  -  DIS  +  XX*(OUT(I)  +  0UT(I+1) ) *.005 
90  DIST( LEVEL)  -  10.*AL0G10( l./DIS) 

RETURN 

END 


QUANTZ 

SUBROUTINE  QOANTZ 

* 

SUBROUTINE  QUANTZ( REQGAN, XBAR, STD, LEVEL , XLEVEL , IOUT, DIST)  1 

C*< 

c* 

*  - 

c* 

SOB  ROUTINE:  QOANTZ 

• 

c* 

PURPOSE: 

c* 

THIS  SUBROUTINE  IS  OSED  TO  SET  OP  THE 

• 

c* 

QUANTIZER  LEVELS  AND  THRESHOLDS 

c* 

n  l 

c* 

EXPLANATION  OF  CALL  VARIABLES: 

c* 

REQGAN-  I  -  REQUIRED  GAIN  OF  THE  QUANTIZER 

c* 

XBAR  -  I  -  MEAN  OF  GUASSIAN  DISTRIBUTION 

*  U. 

c* 

STD  -  I  -  VARIANCE  OF  GAO SIAN  DIST. 

•  - 

c* 

LEVEL  -  0  -  NUMBER  OF  LEVELS  USED  FOR  CODE 

* 

c* 

XLEVEL-  0  -  DISTRIBUTION  THRESHOLDS 

* 

* 

c* 

XOUT  -  0  -  QUANTIZED  OUTPUT  LEVELS 

c* 

DIST  -  I  -  ARRAY  CONTAINING  GAINS  FOR  EACH*  1 

c* 

NUMBER  OF  LEVELS 

* 

c* 

c* 

SUBROUTINES  NEEDED: 

* 

c* 

MDNRIS  -  ROUTINE  TO  FIND  AREA  UNDER  CURVE 

*  i 

c • 

FOR  A  NORMAL ( 0,1)  GAUSSIAN  DIST. 

* 

c* 

(IMSL  LIBRARY) 

* 

c* 

C*< 

i**« 

ft* 

c 

REAL  XDUT<  128), XLEVEL ( 128), DIST(l) 

c 

DETERMINE  THE  NUMBER  OF  LEVELS  REQUIRED  FOR  ERROR  TRANSMISSION  1 

c 

n 

BASED  ON  THE  VALUE  FOR  REQGAN 

< 

5 

II  -  1 

DO  5  1-2,128 

II  -  II  +  1 

IF(RBQGAN.LT.DIST(I) )  GO  TO  6 

CONTINUE 

6 

CONTINUE 

' 

c 

DETERMINE  THE  THRESHOLD  LEVELS 

* 

c 

LEVEL  -  II 

XLEVL  -  LEVEL 

START  -  1. /XLEVL 

NLEVL1  -  LEVEL  -  1 

PERCNT  -  START 

DO  10  1-1 , NLEVL1 

CALL  MDNRIS (PEROrT, XLEVEL ( I ),IER) 

• 

XLEVEL(I)  -  XLEVEL( I) *STD  +  XBAR 

PERCNT  -  PERCNT  +  START 

10 

CONTINUE 

s 

L 

c 

DETERMINE  THE  OUTPUT  LEVELS 

XLEVEL (LEVEL)  -  100000. 

BEGIN  -  START/ 2. 

PERCNT  -  BEGIN 

DO  20  1-1, LEVEL 

CALL  MDNRIS ( PERCNT,  XQUTU ) ,  IER) 

XOUT(I)  -  XOUT(I) *STD  +  XBAR 

PERCNT  -  PERCNT  +  START 

20 

CONTINUE 

RETURN 

END 
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